积性函数求和问题的一种筛法

command_block

2020-01-26 18:43:14

个人记录

stO zzt Orz

本文给出了一种时间,空间复杂度均为 O\Big(\left(\frac{N}{\log N}\right)^{2/3}\Big) 的积性函数求和算法。

前置知识：杜教筛(+贝尔级数+powerful number)

符号约定 :

涉及的积性函数均有 F(1)=1。

对函数 F ,给出范围 N ,对所有 m=\lfloor N/i\rfloor 求解 S_F(m) ,这种操作称之为块筛,有时省略范围。

大 N 一般指题目中给出的 N，而小 n 则用作分析。

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

1.质数的 c 次幂前缀和（优化）

前半部分见 Min_25筛小记。

③

当然，可以考虑利用朴素方法处理一部分 h 的前缀，假设我们预处理 [1,T] ,即 n\leq T 的所有 h(T,\_)。

每次新加质数的时候，会将 T 以内暴力筛除一遍,需要树状数组维护动态前缀和，由于每个数只会被筛除一次，复杂度为 O(T\log T)。

递推时次数为 O\bigg(\sum\limits_{d=1}^{N/T}\frac{\sqrt{N/d}}{logN}\bigg)=O\Big(\frac{N}{logN\sqrt{T}}\Big) ,不过注意到树状数组需要一个 O(\log T) ,这部分复杂度为 O(\frac{N}{\sqrt{T}})。

总的复杂度就是 O(\frac{N}{\sqrt{T}}+TlogT) ,取 T=O((\frac{N}{logN})^{2/3}) 可得复杂度为 O(N^{2/3}log^{1/3}N) 。

④

考虑省去递推运算中的部分树状数组查询操作。

注意到 h(n,k) 中如果 p_k>\sqrt{n} 的话函数值不再改变，可以直接使用静态前缀和。

此时, h(\lfloor n/p_k\rfloor,k-1) 仍在变化的条件是 \sqrt{n/p_k}\geq p_k\rightarrow n\geq p_k^3。

也就是说,在第 k 次转移中，求 h(n,k) 的查询需要树状数组的条件为 p_k^3\leq n。

对于每个 h(n,\_) 有 O(\frac{\sqrt[3]{n}}{\log n}) 个树状数组查询操作,复杂度即为O(\sqrt[3]{n})。

这部分总复杂度是 O\bigg(\sum\limits_{d=1}^{N/T}\sqrt[3]{N/d}\bigg)=O(\frac{N}{T^{2/3}})

这样子的总复杂度是 O\Big(\frac{N}{logN\sqrt{T}}+TlogT\Big) 取 T=O(\frac{N^{2/3}}{log^{4/3}N})可得复杂度为 O(\frac{N^{2/3}}{log^{1/3}N})。

⑤

考虑在 p_k\leq T_2 时暴力递推，这部分复杂度为 O(T_2\sqrt{N}/\log N)。

此时还需筛去的数满足最小素因子大于 T_2 ，且是合数。

笔者在 T\leq 10^8 内验证发现，使用 O(T/\sqrt{T_2}) 拟合效果不错。

参考：当 T_2=50,T=5\times 10^6 时，树状数组修改约为 7\times 10^5 次。（对应 n≈10^{11}）

取 T_2=N^{1/6} ，则复杂度为 O\Big(N^{2/3}/\log N+\frac{N}{logN\sqrt{T}}+T\Big)

令 T=\left(\frac{N}{\log N}\right)^{2/3} 可得复杂度为 O\Big(\left(\frac{N}{\log N}\right)^{2/3}\Big)。

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

2.积性函数求和

问题概述 :

给出一个积性函数 F(n) ,满足如下性质:

给出范围 N ,求出 F 的块筛。

根据 \rm powerful\ number 理论，构造数论函数 G(p)=F(p) ，只需要求出 G 的块筛，就可以在 O(n^{0.6}) 内得到 F 的块筛。

不妨令 G(p^k)=0\ (k>1)。

我们以 k 从大到小的方式,考虑“最小素因子为 p_k 的数”。

设 S_{n,k} 为最小素因子 \geq k ，且素因子次数最高一次的数的集合。

（若素因子次数为 2 ，则 G 函数值为 0 ，可以不考虑）

仿照前面 h 的求解方法,设 r(n,k)=\sum\limits_{i\in S_{n,k}}^nG(i)。

设 D_{n,k} 为最小素因子 =k ，且素因子次数最高一次的数的集合。

按照定义有 S_{n,k}=S_{n,k+1}+D_{n,k}。

可以发现 D_{n,k}=p_kS_{\lfloor n/p_k\rfloor,k+1}。

则 S_{n,k}=S_{n,k+1}+p_kS_{\lfloor n/p_k\rfloor,k+1}

可以据此写出 r 的递推式：

r(n,k)=r(n,k+1)+G(p_k)*\begin{cases}r(\lfloor n/p_k\rfloor,k+1)&(p_k\leq\sqrt{n})\\1&(\sqrt{n}

边界 : s(n,|P|)=1 ,此时只有 G(1)=1 有贡献。

解释 :

当 p_k\leq\sqrt{n} 时,由于不会出现多次质因子，利用积性显然。

当 p_k>\sqrt{n} 时, \lfloor n/p_k\rfloor

当 p_k>n 时, \lfloor n/p_k\rfloor=0 ,无贡献。

现在看来,我们似乎需要枚举所有的质数。

注意到当 p_k>\sqrt{n} 时,只有质数处的 G(p) 会被统计到。

当 p_{t} 为第一个大于 \sqrt{n} 的质数时， s(n,t)=\sum\limits_{\sqrt{n}

现在我们能直接从 $s(n,t)$ 开始往前递推,只需要考虑 $\sqrt{n}$ 以内的质数。

$s$ 的递推式和 $h$ 是相似的,可以使用相同的解法,这里不再赘述了。

不过注意对于 $s(n,k)$ ,当 $p_k>\sqrt{n}$ 的时候,并不是和 $h$ 一样不变。这可以用时间戳解决。

- **例题** : [SP34096 DIVCNTK - Counting Divisors (general)](https://www.luogu.com.cn/problem/SP34096)

这是论文里面的经典题。

为了避免混淆，题目中的 $K$ 使用大写,为了方便我们令 `K++`。

**题意** : 求 $\sum\limits_{i=1}^nd(i^{K-1})$。

在这道题里面,设 $F(n)=d(n^{K-1})$ ,易得这也是个积性函数。

可得 $F(p)=d(p^{K-1})=K;F(p^k)=d(p^{k{K-1}})=(K-1)k+1

直接使用上述算法即可。

这东西常数较小 (远小于洲阁筛,同样精细实现的话比Min_25不相上下)。

在 DIVCNTK 拿到了 LG·Rk1 ,原站 SPOJ·Rk5。

把 k 改成 2,3 均能通过 DIVCNT2 + DIVCNT3 ,而且踩了一堆杜教筛。

目前在 DIVCNT3 甚至比我原来精细实现的杜教筛还快。

#include

#define ull unsigned long long

#define lbit(x) ((x)&-(x))

#define Limit 805000

using namespace std;

int tn,lim,lp,p[Limit>>3],c[Limit],e[Limit];

void getsth()

{

for (int i=2,t;i<=lim;i++){

if (!e[i])c[p[++tn]=i]=1;

for (int j=1;j<=tn&&(t=p[j]*i)<=lim;j++){

e[t]=p[j];

if (i%p[j]==0)break;

if(c[i])c[t]=c[i]+1;

}

}p[tn+1]=Limit;

}

ull N,K,tH[35];

void getPN()

{

tH[0]=1;

for (int k=2;k<=33;k++)

tH[k]=(K-1)*k+1-K*tH[k-1];

}

int tot;

ull t[Limit],x[Limit],tag;

inline void upd(int p,ull x){

while(p<=lim){t[p]+=x;p+=lbit(p);}

}

inline ull qry(int p){

ull sum=0;

while(p){sum+=t[p];p^=lbit(p);}

return sum;

}

void getP()

{

for (tot=1;N>(ull)lim*tot;tot++)x[tot]=N/tot-1;

tot--;

for (int i=1;i<=lim;i++)t[i]=lbit(i);

upd(1,-1);

for (int tp=1;tp<=lp;tp++){

int P=p[tp],u=min(N/((ull)P*P),(ull)tot),

uu=min(u,tot/P);

ull sav=qry(P-1),t=N/P;

for (int i=1;i<=uu;i++)

x[i]-=x[P*i]-sav;

if (t<(1ull<<31))

for (int tt=t,i=uu+1;i<=u;i++)

x[i]-=qry(tt/i)-sav;

else

for (int i=uu+1;i<=u;i++)

x[i]-=qry(t/i)-sav;

if ((ull)P*P<=lim)

for (int j=P*P;j<=lim;j+=P)

if (e[j]==P)upd(j,-1);

}

void sieve()

{

for (lp=1;(ull)p[lp]*p[lp]<=N;lp++);

getP();

memset(t,0,sizeof(ull)*(lim+5));

for (int i=1;i<=tot;i++)x[i]=K*(x[i]-lp);

for (int i=lp+1;p[i]<=lim;i++)upd(p[i],K);

ull tK[35],tag=0;tK[0]=1;

for (int k=1;k<=33;k++)tK[k]=tK[k-1]*K;

for (int i=lp;i;i--){

int P=p[i],u=min(N/((ull)P*P),(ull)tot),

uu=min(u,tot/P);

for (int i=1;i<=uu;i++)

x[i]+=K*(x[P*i]+tag);

ull t=N/P;

if (t<(1ull<<31))

for (int tt=t,i=uu+1;i<=u;i++)

x[i]+=K*qry(tt/i);

else

for (int i=uu+1;i<=u;i++)

x[i]+=K*qry(t/i);

upd(P,K);tag+=K;

if ((ull)P*P<=lim)

for (int j=P*P;j<=lim;j+=P)

if (e[j]==P&&c[j])upd(j,tK[c[j]]);

}for (int i=1;i<=tot;i++)x[i]+=1+tag;

upd(1,1);lim=sqrtl(N);

for (int i=tot+1;i<=lim;i++)x[i]=qry(N/i);

tot=lim;

}

inline ull S(ull d){

if (d<=tot)return x[d];

return qry(N/d);

}

ull ans;

void dfs(ull n,ull num,int t)

{

ans+=num*S(n);

for (int i=t,k;i<=lp;i++){

if ((ull)p[i]*p[i]>N/n)return ;

k=2;

for (ull x=(ull)p[i]*p[i];x*n<=N;k++,x*=p[i])

dfs(n*x,num*tH[k],i+1);

}

#define gL(N) max((int)max(1.75*pow(N/log2(N),0.6666),sqrt(N)+5),100)

void solve()

{

lim=gL(N);if (lim>N)lim=N;

getPN();sieve();

ans=0;dfs(1,1,1);

printf("%llu\n",ans);

}

ull Sn[10050],Sk[10050],mN;

int T;

int main()

{

scanf("%d",&T);

for (int i=1;i<=T;i++){

scanf("%llu%llu",&Sn[i],&Sk[i]);

Sk[i]++;mN=max(mN,Sn[i]);

}lim=gL(mN);if (lim>mN)lim=mN;

getsth();

for (int i=1;i<=T;i++){

N=Sn[i];K=Sk[i];

solve();

}return 0;

}

产品定位方法之产品使用者定位法&使用定位法&分类定位法
电动牙刷如何拆开