首页 > 代码库 > [Spoj]Counting Divisors (cube)
[Spoj]Counting Divisors (cube)
来自FallDream的博客，未经允许，请勿转载，谢谢。
设d(x)表示x的约数个数，求$\sum_{i=1}^{n}d(i^{3})$
There are 5 Input files.
- Input #1: 1≤N≤10000, TL = 1s.
- Input #2: 1≤T≤300, 1≤N≤10^8, TL = 20s.
- Input #3: 1≤T≤75, 1≤N≤10^9, TL = 20s.
- Input #4: 1≤T≤15, 1≤N≤10^10, TL = 20s.
- Input #5: 1≤T≤2, 1≤N≤10^11, TL = 20s.
$i^{3}$的约数个数$d(i^{3})$是一个积性函数，所以转而求$d(x)=\prod{F(pi^{ci})}$,其中$F ( pk ^ {ck} )=3ck+1$
可以直接洲阁筛 学了一天大概懂了 顺便抄了个模板 给个链接吧：大佬博客
-----
gi表示1-i中与前j个质数互质的数字的F之和
fi表示1-i中由小于根号n的后j个质数组成的数字的F之和
容易得出转移方程 $$g[i][j]=g[i][j]-F(pk)g[\frac{i}{pk}][j-1]$$
$$ f[i][j]=f[i][j-1]+\sum_{ck>=1}F(pk^{ck})f[\frac{i}{pk^{ck}}][j]$$
显然i只有根号种取值 对于每个根号n以内的质数都要转移，复杂度$O(\frac{n}{\log n})$
考虑优化，显然$p_{j+1}>i$的时候，g[i][j]=4(3*1+1)
所以当$pj^{2}>i$的时候，g[i][j]=g[i][j-1]+F(pi) 可以不用转移，用的时候补上那一段即可。
之所以把f的状态表示成"后j个",也是出于这个目的
这样的复杂度近似是$O(\frac{n^{\frac{3}{4}}}{logn})$
然后线筛出根号n以内的F[],答案是$f[n]+\sum_{i=1}^{\sqrt{n}}F[i]g[\frac{n}{i}]$
#include<iostream>#include<cstdio>#include<cmath>#define MN 320000#define ll long longusing namespace std;inline ll read(){ ll x = 0; char ch = getchar(); while(ch < ‘0‘ || ch > ‘9‘) ch = getchar(); while(ch >= ‘0‘ && ch <= ‘9‘){x = x * 10 + ch - ‘0‘;ch = getchar();} return x;}int s[MN+5],num=0,last[MN+5],l[MN+5],l0[MN+5],sq,P,N[MN+5];ll f0[MN+5],f[MN+5],g0[MN+5],g[MN+5],d[MN+5],n;bool b[MN+5];void CalcF(){ for(int i=1;i<=sq;++i) f[i]=f0[i]=1; for(int i=P-1;i;--i) { for(int j=1;j<=sq&&l[j]>i;++j) { ll now=(n/j)/s[i]; for(int tms=4;now;now/=s[i],tms+=3) { if(now<=sq) f[j]+=tms*(f0[now]+4*(max(0,N[now]-max(i+1,l0[now])+1))); else f[j]+=tms*(f[n/now]+4*max(0,P-max(i+1,l[n/now]))); } } for(int j=sq;j&&l0[j]>i;--j) { ll now=j/s[i]; for(int tms=4;now;tms+=3,now/=s[i]) f0[j]+=tms*(f0[now]+4*max(0,N[now]-max(i+1,l0[now])+1)); } } for(int i=1;i<=sq;++i) f[i]+=4*(P-l[i]);}void CalcG(){ for(int i=1;i<=sq;++i) g0[i]=i,g[i]=n/i; for(int i=1;i<P;++i) { for(int j=1;j<=sq&&l[j]>i;++j) { ll now=n/j/s[i]; if(now<=sq) g[j]-=g0[now]-max(0,i-l0[now]); else g[j]-=g[n/now]-max(0,i-l[n/now]); } for(int j=sq;j&&l0[j]>i;--j) g0[j]-=g0[j/s[i]]-max(0,i-l0[j/s[i]]); } for(int i=1;i<=sq;++i) g[i]-=P-l[i];}int main(){ d[1]=1; for(int i=2;i<=MN;++i) { if(!b[i]) s[++num]=last[i]=i; for(int j=1;s[j]*i<=MN;++j) { b[s[j]*i]=1,last[s[j]*i]=s[j]; if(i%s[j]==0) break; } int sum=1,tms,p; for(int j=i;j>1;) { tms=0;p=last[j]; for(;j%p==0;j/=p,++tms); sum*=(tms*3+1); } d[i]=sum; N[i]=N[i-1]+(!b[i]); } for(int T=read();T;--T) { n=read();sq=sqrt(n);l[sq+1]=0; for(P=1;1LL*s[P]*s[P]<=n;++P); for(int i=1;i<=sq;++i) for(l0[i]=l0[i-1];1LL*s[l0[i]]*s[l0[i]]<=i;++l0[i]); for(int i=sq;i;--i) for(l[i]=l[i+1];1LL*s[l[i]]*s[l[i]]<=n/i;++l[i]); CalcF();CalcG(); ll ans=f[1]; for(int i=1;i<=sq;++i) ans+=4*d[i]*(g[i]-1); printf("%lld\n",ans); } return 0;}
[Spoj]Counting Divisors (cube)