洲阁筛,一种快速求积性函数前缀和的算法
求$\sum\limits_{i=1}^nF(i)$,其中$F(x)$是积性函数,并且$F(p^c)$是关于$p$的低阶多项式
我们把$1\cdots n$的所有数按是否有$\gt\sqrt n$的质因子分类,那么
$\sum\limits_{i=1}^nF(i)=\sum\limits_{\substack{1\leq i\leq n\\i\,没有\,\gt\sqrt n\,的质因子}}F(i)\left(1+\sum\limits_{\substack{\sqrt n\lt j\leq\left\lfloor\frac ni\right\rfloor\\j\,是质数}}\right)$
当$i\geq\sqrt n$时括号内为$1$,所以我们需要计算两部分
1.对$\forall1\leq i\lt\sqrt n$,计算$\sum\limits_{\substack{\sqrt n\lt j\leq\left\lfloor\frac ni\right\rfloor\\j\,是质数}}F(j)$
2.$\sum\limits_{\substack{\sqrt n\leq i\leq n\\i\,没有\,\gt\sqrt n\,的质因子}}F(i)$
因为$F(p)$是关于$p$的低阶多项式,所以问题转化为求一定范围内的质数幂和
设$g_k(i,j)$表示$[1,j]$中与前$i$个质数互质的数的$k$次幂和,$\leq\sqrt n$的质数有$m$个,我们要求的是$g_k\left(m,\left\lfloor\frac ni\right\rfloor\right)-1$
边界:$g_k(0,i)=\sum\limits_{j=1}^ij^k$,因为$F(p)$是关于$p$的低次多项式,这意味着$k$比较小,那么此式可以快速计算
转移:$g_k(i,j)=g_k(i-1,j)-p_i^kg_k\left(i-1,\left\lfloor\frac j{p_i}\right\rfloor\right)$,减去的是最小质因子$=p_i$的数
因为$j=\left\lfloor\frac n{p_\cdots}\right\rfloor$,所以它的取值只有$O(\sqrt n)$种
当$p_i\gt\left\lfloor\frac j{p_i}\right\rfloor$即$p_i^2\gt j$时,$g_k(i,j)=g_k(i-1,j)-p_i^k$,所以对于一个$j$,我们从小到大枚举$i$,一旦遇到一个$i_0$满足此式,我们以后都不用转移它了,要用到它在$i_1$处的值时减去$\sum\limits_{l=i_0}^{i_1}p_l^k$即可
设$f(i,j)$表示$[1,j]$中仅由$\leq\sqrt n$的后$i$个质数组成的数的$F(x)$之和,我们要求的是$f(m,n)-f(m,\sqrt n-1)$
边界:$f(0,i)=1$
转移:$f(i,j)=f(i-1,j)+\sum\limits_{c\geq1}F(p_i^c)f\left(i-1,\left\lfloor\frac j{p_i^c}\right\rfloor\right)$,加上的是最小质因子$=p_i$的数
同样地,$j$的取值只有$O(\sqrt n)$种
当$p_i^2\gt j$时,$f(i,j)=f(i-1,j)+F(p_i)$,所以我们从大到小枚举质数$p_i$,满足此式时不转移,要用到它的值时加上$[p_i,\min(j,\sqrt n)]$中的质数的$F(p)$即可
总时间复杂度$O\left(\frac{n^{\frac34}}{\log n}\right)$,作者水平不足,证明从略
现在来看spoj上的DIVCNT3,这里$\sigma_0\left(\left(p^c\right)^3\right)=3c+1$,是常数,所以1只用筛质数个数,初始化也很好处理
我的写法可能有一些问题,跑起来挺慢的...
#include<stdio.h> typedef long long ll; const ll T=316500,maxn=635000; ll pr[T+10],md[T+10],c[T+10],d3[T+10],s3[T+10]; bool np[T+10]; void sieve(){ ll i,j,M; M=0; md[1]=1; c[1]=1; d3[1]=1; for(i=2;i<=T;i++){ if(!np[i]){ M++; pr[M]=i; md[i]=i; c[i]=1; d3[i]=4; } for(j=1;j<=M&&i*pr[j]<=T;j++){ np[i*pr[j]]=1; if(i%pr[j]==0){ md[i*pr[j]]=md[i]*pr[j]; c[i*pr[j]]=c[i]+1; d3[i*pr[j]]=d3[i/md[i]]*((c[i]+1)*3+1); break; } md[i*pr[j]]=pr[j]; c[i*pr[j]]=1; d3[i*pr[j]]=d3[i]*4; } } for(i=1;i<=T;i++)s3[i]=s3[i-1]+d3[i]; } const ll mod=1000007; struct map{ ll h[mod],nex[maxn],to[maxn],id[maxn],v[maxn],M; void clear(ll n){ M=0; for(ll i=1;i<=n;i=n/(n/i)+1)h[n/i%mod]=0; } ll&operator[](ll x){ for(ll i=h[x%mod];i;i=nex[i]){ if(id[i]==x)return v[i]; } M++; id[M]=x; nex[M]=h[x%mod]; h[x%mod]=M; return v[M]; } }h; ll bl[maxn],g[maxn],i0[maxn],n,m,C; void getg(){ ll i,j,k; C=0; h.clear(n); for(i=1;i<=n;i=n/(n/i)+1){ C++; bl[C]=n/i; g[C]=n/i; i0[C]=0; h[n/i]=C; } for(i=1;i<=m;i++){ for(j=1;j<=C&&pr[i]*pr[i]<=bl[j];j++){ k=h[bl[j]/pr[i]]; g[j]-=g[k]-(i-1-i0[k]); i0[j]=i; } } for(i=1;i<=C&&pr[m]<=bl[i];i++)g[i]-=m-i0[i]; } ll f[maxn],mp[maxn]; bool us[maxn]; ll ck(ll a,ll b){return a>b?a-b:0;} void getf(){ ll i,j,c,t,p; mp[0]=m; for(i=1;i<=C;i++){ us[i]=0; f[i]=1; for(j=mp[i-1];pr[j]>bl[i];j--); mp[i]=j; } for(i=m;i>0;i--){ for(j=1;j<=C&&pr[i]*pr[i]<=bl[j];j++){ if(!us[j]){ f[j]+=ck(mp[j],i)*4; us[j]=1; } for(c=1,t=pr[i];t<=bl[j];c++,t*=pr[i]){ p=h[bl[j]/t]; f[j]+=(f[p]+(us[p]?0:ck(mp[p],i)*4))*(3*c+1); } } } for(i=1;i<=C;i++){ if(!us[i])f[i]+=mp[i]*4; } } void work(){ ll i,sq,ans; scanf("%lld",&n); if(n<=T){ printf("%lld\n",s3[n]); return; } for(m=1;pr[m]*pr[m]<=n;m++); m--; getg(); getf(); for(sq=1;sq*sq<n;sq++); sq--; ans=f[1]-f[h[sq]]; for(i=1;i<=sq;i++)ans+=d3[i]*((g[h[n/i]]-1)*4+1); printf("%lld\n",ans); } int main(){ sieve(); ll cas; scanf("%lld",&cas); while(cas--)work(); }
DIVCNTK被卡常了,鈤