简介:
对于积性函数f(x)f(x)f(x)
求∑i=1i≤nf(i)(其中n≤1011左右)\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)=A0+A1P+A2P2+……F(P)=A_0+A_1P+A_2P^2+……F(P)=A0+A1P+A2P2+……。
并且项数不能太多(因为每一项都要手写公式求和)
算法概述:
在简介中已经写到,对于每一项都要分别求和,因此假设我们现在求的是第k+1项,即f(P)=AkPkf(P)=A_kP^kf(P)=AkPk
设g(i,j)=∑k为质数,或k的最小质因子>Pjf(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)} 22i∗(i+1)……
然后需要转移:
首先,要求i≥prj∗prji\geq pr_j*pr_ji≥prj∗prj
(否则g(i,j)=g(i,j−1)g(i,j)=g(i,j-1)g(i,j)=g(i,j−1))
当然,暴力算还是会T的。
但是由于一个性质:(nab\frac{\frac {n} {a}} {b}ban的个数=nc\frac {n} {c}cn的个数,即不超过2n2\sqrt n2n种)
就可以预处理出可能的iPj\frac {i} {P_j}Pji的值。然后就能快速算了。
然后,设
(当i<prji< pr_ji<prj时,S(i,j)=0S(i,j)=0S(i,j)=0,因此,下面默认i≥prji\geq pr_ji≥prj)
显然,其中的质数部分为:g(i,∣P∣)−∑k=1j−1f(Pk)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);
}