简介:
对于积性函数 f ( x ) f(x) f(x)
求 ∑ i = 1 i ≤ n f ( i ) ( 其 中 n ≤ 1 0 11 左 右 ) \sum_{i=1}^{i\leq n} f(i)(其中n\leq 10^{11}左右) i=1∑i≤nf(i)(其中n≤1011左右)
必须满足的条件是:
当p为质数时,
F
(
P
)
F(P)
F(P)必须能表示为一个多项式的形式,即
F
(
P
)
=
A
0
+
A
1
P
+
A
2
P
2
+
…
…
F(P)=A_0+A_1P+A_2P^2+……
F(P)=A0+A1P+A2P2+……。
并且项数不能太多(因为每一项都要手写公式求和)
算法概述:
在简介中已经写到,对于每一项都要分别求和,因此假设我们现在求的是第k+1项,即 f ( P ) = A k P k f(P)=A_kP^k f(P)=AkPk
设 g ( i , j ) = ∑ k 为 质 数 , 或 k 的 最 小 质 因 子 > P j f ( k ) g(i,j)=\sum_{k为质数,或k的最小质因子>P_j}f(k) g(i,j)=k为质数,或k的最小质因子>Pj∑f(k)
那么,显然
g
(
i
,
0
)
g(i,0)
g(i,0)就是所有数的k次方和。
这玩意就得考手推公式了:比如0次方和=i,1次方和=
i
∗
(
i
+
1
)
2
\frac {i*(i+1)} 2
2i∗(i+1)……
然后需要转移:
首先,要求
i
≥
p
r
j
∗
p
r
j
i\geq pr_j*pr_j
i≥prj∗prj
(否则
g
(
i
,
j
)
=
g
(
i
,
j
−
1
)
g(i,j)=g(i,j-1)
g(i,j)=g(i,j−1))
当然,暴力算还是会T的。
但是由于一个性质:(
n
a
b
\frac{\frac {n} {a}} {b}
ban的个数=
n
c
\frac {n} {c}
cn的个数,即不超过
2
n
2\sqrt n
2n种)
就可以预处理出可能的
i
P
j
\frac {i} {P_j}
Pji的值。然后就能快速算了。
然后,设
(当
i
<
p
r
j
i< pr_j
i<prj时,
S
(
i
,
j
)
=
0
S(i,j)=0
S(i,j)=0,因此,下面默认
i
≥
p
r
j
i\geq pr_j
i≥prj)
显然,其中的质数部分为:
g
(
i
,
∣
P
∣
)
−
∑
k
=
1
j
−
1
f
(
P
k
)
g(i,|P|)-\sum_{k=1}^{j-1}f(P_k)
g(i,∣P∣)−∑k=1j−1f(Pk)
然后也可以转移:
套路
1、对于一些题目,可能在某些特定位置不满足条件,此时可以先忽略,最后再特判这些位置。
2、有些函数在多项式中系数不为1,这时最好先忽略系数,最后来乘。
板子题
LOJ6053简单的函数
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#define SF scanf
#define PF printf
#define MAXN 200010
#define MOD 1000000007
using namespace std;
typedef long long ll;
ll sqr,n;
bool isprime[MAXN];
int primes[MAXN],tot;
ll sum[MAXN];
void prepare(){
for(int i=2;i<=sqr;i++){
if(isprime[i]==0){
primes[++tot]=i;
sum[tot]=(sum[tot-1]+i)%MOD;
}
for(int j=1;1ll*i*primes[j]<=sqr;j++){
isprime[i*primes[j]]=1;
if(i%primes[j]==0)
break;
}
}
}
ll w[MAXN],id1[MAXN],id2[MAXN],g[MAXN],h[MAXN];
int cnt;
void get_g(){
ll las;
for(ll i=1;i<=n;i=las+1){
las=n/(n/i);
w[++cnt]=n/i;
if(w[cnt]<=sqr)
id1[w[cnt]]=cnt;
else
id2[n/w[cnt]]=cnt;
g[cnt]=(w[cnt]%MOD+2)*(w[cnt]%MOD-1)%MOD*((MOD+1ll)/2ll)%MOD;
h[cnt]=w[cnt]%MOD-1;
}
for(int j=1;j<=tot;j++)
for(int i=1;i<=cnt&&1ll*primes[j]*primes[j]<=w[i];i++){
ll k=w[i]/primes[j];
if(k<=sqr)
k=id1[k];
else
k=id2[n/k];
g[i]=(g[i]-1ll*primes[j]*(g[k]-sum[j-1])%MOD+MOD)%MOD;
h[i]=(h[i]-1ll*(h[k]-(j-1))%MOD+MOD)%MOD;
}
}
int S(ll x,int j){
if(x<=1||primes[j]>x)
return 0;
int k1;
if(x<=sqr)
k1=id1[x];
else
k1=id2[n/x];
ll res=((g[k1]-h[k1]-sum[j-1]+j-1)%MOD+MOD)%MOD;
if(j==1)
res=(res+2)%MOD;
for(int k=j;k<=tot&&1ll*primes[k]*primes[k]<=x;k++){
ll t=primes[k];
for(int e=1;1ll*primes[k]*t<=x;e++,t*=primes[k]){
res=((res+1ll*(primes[k]^e)*S(x/t,k+1)%MOD)%MOD+(primes[k]^(e+1)))%MOD;
}
}
return res;
}
int main(){
SF("%lld",&n);
sqr=sqrt(n);
prepare();
// PF("[%lld %d]",sqr,tot);
get_g();
PF("%d",(S(n,1)+1)%MOD);
}