[SPOJ]DIVCNT3

别人写的讲得挺好的博客

洲阁筛，一种快速求积性函数前缀和的算法

求$sumlimits_{i=1}^nF(i)$，其中$F(x)$是积性函数，并且$F(p^c)$是关于$p$的低阶多项式

我们把$1cdots n$的所有数按是否有$gtsqrt n$的质因子分类，那么

$sumlimits_{i=1}^nF(i)=sumlimits_{substack{1leq ileq n\i,没有,gtsqrt n,的质因子}}F(i)left(1+sumlimits_{substack{sqrt nlt jleqleftlfloorfrac nirightrfloor\j,是质数}}right)$

当$igeqsqrt n$时括号内为$1$，所以我们需要计算两部分

1.对$forall1leq iltsqrt n$，计算$sumlimits_{substack{sqrt nlt jleqleftlfloorfrac nirightrfloor\j,是质数}}F(j)$

2.$sumlimits_{substack{sqrt nleq ileq n\i,没有,gtsqrt n,的质因子}}F(i)$

因为$F(p)$是关于$p$的低阶多项式，所以问题转化为求一定范围内的质数幂和

设$g_k(i,j)$表示$[1,j]$中与前$i$个质数互质的数的$k$次幂和，$leqsqrt n$的质数有$m$个，我们要求的是$g_kleft(m,leftlfloorfrac nirightrfloorright)-1$

边界：$g_k(0,i)=sumlimits_{j=1}^ij^k$，因为$F(p)$是关于$p$的低次多项式，这意味着$k$比较小，那么此式可以快速计算

转移：$g_k(i,j)=g_k(i-1,j)-p_i^kg_kleft(i-1,leftlfloorfrac j{p_i}rightrfloorright)$，减去的是最小质因子$=p_i$的数

因为$j=leftlfloorfrac n{p_cdots}rightrfloor$，所以它的取值只有$O(sqrt n)$种

当$p_igtleftlfloorfrac j{p_i}rightrfloor$即$p_i^2gt j$时，$g_k(i,j)=g_k(i-1,j)-p_i^k$，所以对于一个$j$，我们从小到大枚举$i$，一旦遇到一个$i_0$满足此式，我们以后都不用转移它了，要用到它在$i_1$处的值时减去$sumlimits_{l=i_0}^{i_1}p_l^k$即可

设$f(i,j)$表示$[1,j]$中仅由$leqsqrt n$的后$i$个质数组成的数的$F(x)$之和，我们要求的是$f(m,n)-f(m,sqrt n-1)$

边界：$f(0,i)=1$

转移：$f(i,j)=f(i-1,j)+sumlimits_{cgeq1}F(p_i^c)fleft(i-1,leftlfloorfrac j{p_i^c}rightrfloorright)$，加上的是最小质因子$=p_i$的数

同样地，$j$的取值只有$O(sqrt n)$种

当$p_i^2gt j$时，$f(i,j)=f(i-1,j)+F(p_i)$，所以我们从大到小枚举质数$p_i$，满足此式时不转移，要用到它的值时加上$[p_i,min(j,sqrt n)]$中的质数的$F(p)$即可

总时间复杂度$Oleft(frac{n^{frac34}}{log n}right)$，作者水平不足，证明从略

现在来看spoj上的DIVCNT3，这里$sigma_0left(left(p^cright)^3right)=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("%lldn",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("%lldn",ans);
}
int main(){
sieve();
ll cas;
scanf("%lld",&cas);
while(cas--)work();
}

DIVCNTK被卡常了，鈤