[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(); }

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#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被卡常了，鈤