算法总结积性函数相关
Posted zsnuo
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了算法总结积性函数相关相关的知识，希望对你有一定的参考价值。
【线性筛】
〖模板代码〗
[线性筛质数]
1 int main() 2 { 3 n=read(); 4 for(int i=2;i<=n;i++) 5 { 6 if(!f[i])pri[++cnt]=i; 7 for(int j=1;j<=cnt;j++) 8 { 9 if(i*pri[j]>n)break; 10 f[i*pri[j]]=1; 11 if(i%pri[j]==0)break; 12 } 13 } 14 }
[线性求约数和]
1 void pre() 2 { 3 miu[1]=1;f[1].first=1; 4 for(int i=2;i<=mx;i++) 5 { 6 if(!np[i]){pri[++cnt]=i;d[i]=i;f[i].first=i+1;miu[i]=-1;} 7 for(int j=1;i*pri[j]<=mx;j++) 8 { 9 np[i*pri[j]]=true; 10 if(i%pri[j]==0) 11 { 12 miu[i*pri[j]]=0; 13 d[i*pri[j]]=d[i]*pri[j]; 14 f[i*pri[j]].first=f[i].first+f[i/d[i]].first*d[i*pri[j]]; 15 break; 16 } 17 miu[i*pri[j]]=-miu[i]; 18 d[i*pri[j]]=pri[j]; 19 f[i*pri[j]].first=f[i].first*f[pri[j]].first; 20 } 21 } 22 for(int i=1;i<=mx;i++)f[i].second=i; 23 }
〖相关题目〗
1.【Technocup 2018 - Elimination Round 2】F. Paths
题意：给定数字n，建立一个无向图。对于所有1~n之间的数字，当数字gcd(u,v)≠1时将u、v连一条边，边权为1。d(u, v)表示u到v的最短路，求所有d(u, v)的和，其中1 ≤ u < v ≤ n。
分析：Zsnuoの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #define LL long long 5 using namespace std; 6 const int N=1e7+5; 7 int n,m,tot,now,pri[N],p[N],phi[N],num[N],sum[N]; 8 LL one,two,three; 9 int main() 10 { 11 scanf("%d",&n); 12 phi[1]=1; 13 for(int i=2;i<=n;i++) 14 { 15 if(!p[i]){p[i]=pri[++tot]=i;phi[i]=i-1;} 16 for(int j=1;j<=tot;j++) 17 { 18 if(i*pri[j]>n)break; 19 p[i*pri[j]]=pri[j]; 20 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 21 else phi[i*pri[j]]=phi[i]*(pri[j]-1); 22 } 23 } 24 for(int i=2;i<=n;i++)one+=i-1-phi[i]; 25 for(int i=2;i<=n;i++)num[p[i]]++; 26 for(int i=2;i<=n;i++)sum[i]=sum[i-1]+num[i]; 27 for(int i=2;i<=n;i++)two+=1ll*num[i]*sum[n/i]; 28 for(int i=2;i<=n;i++)if(1ll*p[i]*p[i]<=n)two--; 29 two=two/2-one;m=n-1; 30 for(int i=tot;i>=1;i--) 31 if(pri[i]*2>n)m--; 32 else break; 33 three=1ll*m*(m-1)/2-one-two; 34 printf("%I64d\\n",one+two*2+three*3); 35 return 0; 36 }
2.【51nod1584】加权约数和
题意：见原题
分析：KsClaの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=1e6; 8 const int mod=1e9+7; 9 int T,n,tot,now,pri[N+5],miu[N+5]; 10 int mn[N+5],h[N+5],p[N+5],g[N+5],f[N+5]; 11 bool np[N+5]; 12 int read() 13 { 14 int x=0,f=1;char c=getchar(); 15 while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1;c=getchar();} 16 while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘;c=getchar();} 17 return x*f; 18 } 19 int main() 20 { 21 h[1]=p[1]=miu[1]=g[1]=1; 22 for(int i=2;i<=N;i++) 23 { 24 if(!np[i]) 25 { 26 pri[++tot]=i; 27 miu[i]=-1;mn[i]=i;h[i]=i+1; 28 p[i]=(1ll*i*i%mod+i+1)%mod; 29 } 30 for(int j=1;i*pri[j]<=N;j++) 31 { 32 now=i*pri[j];np[now]=true; 33 if(i%pri[j]==0) 34 { 35 miu[now]=0;mn[now]=mn[i]*pri[j]; 36 h[now]=(h[i]+1ll*h[i/mn[i]]*mn[now]%mod)%mod; 37 p[now]=(1ll*(mn[now]+mn[i])%mod*mn[now]%mod*p[i/mn[i]]%mod+p[i])%mod; 38 break; 39 } 40 miu[now]=-miu[i];mn[now]=pri[j]; 41 h[now]=1ll*h[i]*h[pri[j]]%mod; 42 p[now]=1ll*p[i]*p[pri[j]]%mod; 43 } 44 } 45 for(int i=2;i<=N;i++)p[i]=1ll*p[i]*i%mod; 46 for(int i=2;i<=N;i++)p[i]=(p[i]+p[i-1])%mod; 47 for(int i=2;i<=N;i++)miu[i]=1ll*(miu[i]+mod)%mod*i%mod*i%mod; 48 for(int i=2;i<=N;i++)g[i]=(g[i-1]+h[i])%mod; 49 for(int i=2;i<=N;i++)g[i]=1ll*i*g[i]%mod*h[i]%mod; 50 for(int i=1;i<=N;i++) 51 for(int j=1;i*j<=N;j++) 52 f[i*j]=(f[i*j]+1ll*miu[i]*g[j]%mod)%mod; 53 for(int i=2;i<=N;i++)f[i]=(f[i-1]+f[i])%mod; 54 T=read(); 55 for(int i=1;i<=T;i++) 56 { 57 n=read(); 58 n=(2*f[n]%mod-p[n]+mod)%mod; 59 printf("Case #%d: %d\\n",i,n); 60 } 61 return 0; 62 }
【杜教筛】
〖相关资料〗
《浅谈一类积性函数的前缀和》 《论逗逼的自我修养之寒假颓废记》
〖相关题目〗
题意：求φ的前缀和。
分析：ONION_CYCの博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=5e6; 8 const int mod=1e9+7; 9 const int inv=(mod+1)/2; 10 int s,tot,pri[N+5],phi[N+5],a[N+5],b[N+5]; 11 bool isp[N+5]; 12 LL n; 13 int memory(LL x) 14 { 15 if(x<=s)return a[x]; 16 else return b[n/x]; 17 } 18 void add(LL x,int ans) 19 { 20 if(x<=s)a[x]=ans; 21 else b[n/x]=ans; 22 } 23 int solve(LL n) 24 { 25 if(n<=N)return phi[n]; 26 if(~(memory(n)))return memory(n); 27 int ans=1ll*(n%mod)*((n+1)%mod)%mod*inv%mod; 28 LL pos; 29 for(LL i=2;i<=n;i=pos+1) 30 { 31 pos=n/(n/i); 32 ans=(ans-1ll*(pos-i+1)%mod*solve(n/i)%mod+mod)%mod; 33 } 34 add(n,ans); 35 return ans; 36 } 37 int main() 38 { 39 memset(a,-1,sizeof(a)); 40 memset(b,-1,sizeof(b)); 41 scanf("%lld",&n); 42 s=sqrt(n);phi[1]=1; 43 for(int i=2;i<=N;i++) 44 { 45 if(!isp[i]){pri[++tot]=i;phi[i]=i-1;} 46 for(int j=1;j<=tot;j++) 47 { 48 if(i*pri[j]>N)break; 49 isp[i*pri[j]]=true; 50 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 51 else phi[i*pri[j]]=phi[i]*(pri[j]-1); 52 } 53 phi[i]=(phi[i]+phi[i-1])%mod; 54 } 55 if(n<=N)printf("%d",phi[n]); 56 else printf("%d",solve(n)); 57 return 0; 58 }
2.【51nod1244】莫比乌斯函数之和
题意：求μ的前缀和
分析：无
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=5e6; 8 int s,tot,pri[N+5],miu[N+5],a[N+5]; 9 bool np[N+5]; 10 LL nn,L,R; 11 LL solve(LL n) 12 { 13 if(n<=N)return miu[n]; 14 if(~a[nn/n])return a[nn/n]; 15 LL ans=1,pos; 16 for(LL i=2;i<=n;i=pos+1) 17 { 18 pos=n/(n/i); 19 ans=ans-(pos-i+1)*solve(n/i); 20 } 21 return a[nn/n]=ans; 22 } 23 int main() 24 { 25 miu[1]=1; 26 for(int i=2;i<=N;i++) 27 { 28 if(!np[i]){miu[i]=-1;pri[++tot]=i;} 29 for(int j=1;i*pri[j]<=N;j++) 30 { 31 np[i*pri[j]]=true; 32 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 33 miu[i*pri[j]]=-miu[i]; 34 } 35 miu[i]+=miu[i-1]; 36 } 37 scanf("%lld%lld",&L,&R); 38 nn=L-1;memset(a,-1,sizeof(a));L=solve(L-1); 39 nn=R;memset(a,-1,sizeof(a));R=solve(R); 40 printf("%lld",R-L); 41 return 0; 42 }
3.【bzoj3944】Sum
题意：求φ和μ的前缀和
分析：无
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=2e6; 8 LL T,n,nn,tot,pri[N+5]; 9 LL miu[N+5],phi[N+5],mmiu[N+5],mphi[N+5]; 10 bool np[N+5]; 11 struct node{LL phi,miu;}ans,tmp; 12 node solve(LL n) 13 { 14 if(n<=N)return (node){phi[n],miu[n]}; 15 if(~mphi[nn/n])return (node){mphi[nn/n],mmiu[nn/n]}; 16 LL ansp=1ll*n*(n+1)/2,ansm=1,pos; 17 for(LL i=2;i<=n;i=pos+1) 18 { 19 pos=n/(n/i);tmp=solve(n/i); 20 ansp=ansp-(pos-i+1)*tmp.phi; 21 ansm=ansm-(pos-i+1)*tmp.miu; 22 } 23 mphi[nn/n]=ansp;mmiu[nn/n]=ansm; 24 return (node){ansp,ansm}; 25 } 26 void work() 27 { 28 scanf("%lld",&n);nn=n; 29 memset(mphi,-1,sizeof(mphi)); 30 ans=solve(n); 31 printf("%lld %lld\\n",ans.phi,ans.miu); 32 } 33 int main() 34 { 35 phi[1]=miu[1]=1; 36 for(int i=2;i<=N;i++) 37 { 38 if(!np[i]){pri[++tot]=i;miu[i]=-1;phi[i]=i-1;} 39 for(int j=1;i*pri[j]<=N;j++) 40 { 41 np[i*pri[j]]=true; 42 if(i%pri[j]==0) 43 { 44 miu[i*pri[j]]=0; 45 phi[i*pri[j]]=phi[i]*pri[j]; 46 break; 47 } 48 miu[i*pri[j]]=-miu[i]; 49 phi[i*pri[j]]=phi[i]*(pri[j]-1); 50 } 51 miu[i]+=miu[i-1]; 52 phi[i]+=phi[i-1]; 53 } 54 scanf("%lld",&T); 55 while(T--)work(); 56 return 0; 57 }
4.【51nod1238】最小公倍数之和V3
题意：输出小于等于N的所有数，两两之间的最小公倍数之和。
分析：jiry_2の博客
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=5e6; 8 const int mod=1e9+7; 9 const int inv2=(mod+1)/2; 10 const int inv6=(mod+1)/6; 11 int tot,now,pri[N+5],mn[N+5]; 12 bool np[N+5]; 13 LL n,nn,sum,phi[N+5],a[N+5]; 14 LL calc1(LL a) 15 { 16 LL suma=a%mod; 17 return suma*(suma+1)%mod*inv2%mod; 18 } 19 LL calc2(LL a,LL b) 20 { 21 LL suma=(a-1)%mod; 22 suma=suma*(suma+1)%mod*((2*suma+1)%mod)%mod*inv6%mod; 23 LL sumb=b%mod; 24 sumb=sumb*(sumb+1)%mod*((2*sumb+1)%mod)%mod*inv6%mod; 25 return (sumb-suma+mod)%mod; 26 } 27 LL calc3(LL a,LL b) 28 { 29 LL suma=(a-1)%mod; 30 suma=suma*(suma+1)%mod*inv2%mod; 31 suma=suma*suma%mod; 32 LL sumb=b%mod; 33 sumb=sumb*(sumb+1)%mod*inv2%mod; 34 sumb=sumb*sumb%mod; 35 return (sumb-suma+mod)%mod; 36 } 37 LL solve(LL n) 38 { 39 if(n<=N)return phi[n]; 40 if(~a[nn/n])return a[nn/n]; 41 LL ans=0,pos; 42 for(LL i=1;i<=n;i=pos+1) 43 { 44 pos=n/(n/i); 45 ans=(ans+calc1(n/i)*calc3(i,pos))%mod; 46 } 47 for(LL i=2;i<=n;i=pos+1) 48 { 49 pos=n/(n/i); 50 ans=(ans-solve(n/i)*calc2(i,pos)%mod+mod)%mod; 51 } 52 return a[nn/n]=ans; 53 } 54 int main() 55 { 56 scanf("%lld",&n); 57 nn=n;phi[1]=1; 58 for(int i=2;i<=N;i++) 59 { 60 if(!np[i]) 61 { 62 pri[++tot]=i; 63 phi[i]=(1ll*i*(i-1)%mod+1)%mod; 64 mn[i]=i; 65 } 66 for(int j=1;i*pri[j]<=N;j++) 67 { 68 now=i*pri[j]; 69 np[now]=true; 70 if(i%pri[j]==0) 71 { 72 mn[now]=mn[i]*pri[j]; 73 phi[now]=(phi[i]+phi[i/mn[i]]*mn[now]%mod*(mn[i]*(pri[j]-1)%mod))%mod; 74 break; 75 } 76 phi[now]=phi[i]*(1ll*pri[j]*(pri[j]-1)%mod+1)%mod; 77 mn[now]=pri[j]; 78 } 79 } 80 for(int i=1;i<=N;i++)phi[i]=phi[i]*i%mod; 81 for(int i=2;i<=N;i++)phi[i]=(phi[i]+phi[i-1])%mod; 82 memset(a,-1,sizeof(a)); 83 printf("%lld\\n",solve(n)); 84 return 0; 85 }
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=5e6; 8 const int mod=1e9+7; 9 const int inv2=(mod+1)/2; 10 const int inv6=(mod+1)/6; 11 int tot,pri[N+5]; 12 LL n,nn; 13 LL phi[N+5],a[N+5]; 14 bool np[N+5]; 15 LL calc(LL a,LL b) 16 { 17 LL suma=(a-1)%mod; 18 suma=suma*(suma+1)%mod*((2*suma+1)%mod)%mod*inv6%mod; 19 LL sumb=b%mod; 20 sumb=sumb*(sumb+1)%mod*((2*sumb+1)%mod)%mod*inv6%mod; 21 return (sumb-suma+mod)%mod; 22 } 23 LL solve(LL n) 24 { 25 if(n<=N)return phi[n]; 26 if(~a[nn/n])return a[nn/n]; 27 LL ans=(n%mod)*((n+1)%mod)%mod*inv2%mod,pos; 28 ans=ans*ans%mod; 29 for(LL i=2;i<=n;i=pos+1) 30 { 31 pos=n/(n/i); 32 ans=(ans-solve(n/i)*calc(i,pos)%mod+mod)%mod; 33 } 34 return a[nn/n]=ans; 35 } 36 int main() 37 { 38 scanf("%lld",&n); 39 nn=n;phi[1]=1; 40 for(int i=2;i<=N;i++) 41 { 42 if(!np[i]){pri[++tot]=i;phi[i]=i-1;} 43 for(int j=1;i*pri[j]<=N;j++) 44 { 45 np[i*pri[j]]=true; 46 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 47 phi[i*pri[j]]=phi[i]*(pri[j]-1); 48 } 49 phi[i]=phi[i]*i%mod*i%mod; 50 phi[i]=(phi[i]+phi[i-1])%mod; 51 } 52 memset(a,-1,sizeof(a)); 53 LL pos,ans=0,sum; 54 for(LL i=1;i<=n;i=pos+1) 55 { 56 pos=n/(n/i); 57 sum=(solve(pos)-solve(i-1)+mod)%mod; 58 // printf(" %lld %lld\\n",pos,solve(pos)); 59 sum=sum*(((n/i)%mod)*((n/i+1)%mod)%mod*inv2%mod)%mod; 60 ans=(ans+sum)%mod; 61 } 62 printf("%lld\\n",ans); 63 // for(int i=1;i<=n;i++)printf("%lld\\n",phi[i]); 64 return 0; 65 }
5.【51nod1237】最大公约数之和
题意：输出小于等于N的所有数，两两之间的最大公约数之和。
分析：无
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=5e6; 8 const int mod=1e9+7; 9 const int inv=(mod+1)/2; 10 int tot,pri[N+5]; 11 LL n,nn,ans,pos,phi[N+5],a[N+5]; 12 bool np[N+5]; 13 LL calc(LL a,LL b) 14 { 15 LL suma=(a-1)%mod; 16 suma=suma*(suma+1)%mod*inv%mod; 17 LL sumb=b%mod; 18 sumb=sumb*(sumb+1)%mod*inv%mod; 19 return (sumb-suma+mod)%mod; 20 } 21 LL solve(LL n) 22 { 23 if(n<=N)return phi[n]; 24 if(~a[nn/n])return a[nn/n]; 25 LL ans=(n%mod)*((n+1)%mod)%mod*inv%mod,pos; 26 for(LL i=2;i<=n;i=pos+1) 27 { 28 pos=n/(n/i); 29 ans=(ans-(pos-i+1)%mod*solve(n/i)%mod+mod)%mod; 30 } 31 return a[nn/n]=ans; 32 } 33 int main() 34 { 35 phi[1]=1; 36 for(int i=2;i<=N;i++) 37 { 38 if(!np[i]){pri[++tot]=i;phi[i]=i-1;} 39 for(int j=1;i*pri[j]<=N;j++) 40 { 41 np[i*pri[j]]=true; 42 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 43 phi[i*pri[j]]=phi[i]*(pri[j]-1); 44 } 45 phi[i]=(phi[i]+phi[i-1])%mod; 46 } 47 scanf("%lld",&n);nn=n; 48 memset(a,-1,sizeof(a)); 49 for(LL i=1;i<=n;i=pos+1) 50 { 51 pos=n/(n/i); 52 ans=(ans+calc(i,pos)*solve(n/i)%mod)%mod; 53 } 54 n%=mod;n=n*(n+1)%mod*inv%mod; 55 ans=(ans*2%mod-n+mod)%mod; 56 printf("%lld",ans); 57 return 0; 58 }
6.【51nod1227】平均最小公倍数
题意：见原题。
分析：无
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=1e6; 8 const int mod=1e9+7; 9 const int inv2=(mod+1)/2; 10 const int inv6=(mod+1)/6; 11 int tot,pri[N+5]; 12 LL a,b,nn,pos,sum,phi[N+5],m[N+5]; 13 bool np[N+5]; 14 LL calc(LL a,LL b) 15 { 16 LL suma=(a-1)%mod; 17 suma=suma*(suma+1)%mod*inv2%mod; 18 LL sumb=b%mod; 19 sumb=sumb*(sumb+1)%mod*inv2%mod; 20 return (sumb-suma+mod)%mod; 21 } 22 LL solve(LL n) 23 { 24 if(n<=N)return phi[n]; 25 if(~m[nn/n])return m[nn/n]; 26 LL ans=n%mod,pos; 27 ans=ans*(ans+1)%mod*((ans*2+1)%mod)%mod*inv6%mod; 28 for(LL i=2;i<=n;i=pos+1) 29 { 30 pos=n/(n/i); 31 sum=calc(i,pos)*solve(n/i)%mod; 32 ans=(ans-sum+mod)%mod; 33 } 34 return m[nn/n]=ans; 35 } 36 LL work(LL n) 37 { 38 LL ans=0; 39 for(LL i=1;i<=n;i=pos+1) 40 { 41 pos=n/(n/i); 42 sum=((pos-i+1)%mod)*solve(n/i)%mod; 43 ans=(ans+sum)%mod; 44 // printf("%lld\\n",ans); 45 } 46 return (ans+n)%mod*inv2%mod; 47 } 48 int main() 49 { 50 phi[1]=1; 51 for(int i=2;i<=N;i++) 52 { 53 if(!np[i]){pri[++tot]=i;phi[i]=i-1;} 54 for(int j=1;i*pri[j]<=N;j++) 55 { 56 np[i*pri[j]]=true; 57 if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;} 58 phi[i*pri[j]]=phi[i]*(pri[j]-1); 59 } 60 phi[i]=phi[i]*i%mod; 61 phi[i]=(phi[i]+phi[i-1])%mod; 62 } 63 scanf("%lld%lld",&a,&b); 64 memset(m,-1,sizeof(m));nn=a;a=work(a-1); 65 memset(m,-1,sizeof(m));nn=b;b=work(b); 66 printf("%lld",(b-a+mod)%mod); 67 return 0; 68 }
7.【bzoj4176】Lucas的数论
题意：见原题。
分析：PoPoQQQの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=6e6; 8 const int mod=1e9+7; 9 int S,tot,pri[N+5],miu[N+5],a[N+5]; 10 int n,nn,pos,n2,pos2,sum,ans; 11 bool np[N+5]; 12 int findmiu(int n) 13 { 14 if(n<=S)return miu[n]; 15 if(~a[nn/n])return a[nn/n]; 16 int ans=1,pos; 17 for(int i=2;i<=n;i=pos+1) 18 { 19 pos=n/(n/i); 20 ans=(ans-(LL)(pos-i+1)*findmiu(n/i)%mod+mod)%mod; 21 } 22 return a[nn/n]=ans; 23 } 24 int solve(int a,int b) 25 { 26 a=findmiu(a-1);b=findmiu(b); 27 return (b-a+mod)%mod; 28 } 29 int main() 30 { 31 scanf("%d",&n); 32 S=ceil(pow(n,0.75)-1e-7)+1e-7; 33 miu[1]=1; 34 for(int i=2;i<=S;i++) 35 { 36 if(!np[i]){pri[++tot]=i;miu[i]=-1;} 37 for(int j=1;i*pri[j]<=S;j++) 38 { 39 np[i*pri[j]]=true; 40 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 41 miu[i*pri[j]]=-miu[i]; 42 } 43 miu[i]+=miu[i-1]; 44 } 45 nn=n;memset(a,-1,sizeof(a)); 46 for(int i=1;i<=n;i=pos+1) 47 { 48 pos=n/(n/i);n2=n/i;sum=0; 49 for(int j=1;j<=n2;j=pos2+1) 50 { 51 pos2=n2/(n2/j); 52 sum=(sum+(LL)(pos2-j+1)*(n2/j)%mod)%mod; 53 } 54 sum=(LL)sum*sum%mod; 55 ans=(ans+(LL)solve(i,pos)*sum%mod)%mod; 56 } 57 printf("%d",ans); 58 return 0; 59 }
8.【51nod1220】约数之和
题意：见原题。
分析：无。
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=6e6; 8 const int mod=1e9+7; 9 int n,S,nn,tot,pos,ans; 10 int pri[N+5],miu[N+5],a[N+5]; 11 bool np[N+5]; 12 void pre() 13 { 14 S=ceil(pow(n,0.75)-1e-7)+1e-7; 15 miu[1]=1; 16 for(int i=2;i<=S;i++) 17 { 18 if(!np[i]){pri[++tot]=i;miu[i]=-1;} 19 for(int j=1;i*pri[j]<=S;j++) 20 { 21 np[i*pri[j]]=true; 22 if(i%pri[j]==0){miu[i*pri[j]]=0;break;} 23 miu[i*pri[j]]=-miu[i]; 24 } 25 miu[i]=1ll*miu[i]*i%mod; 26 miu[i]=(miu[i-1]+miu[i])%mod; 27 miu[i]=(miu[i]+mod)%mod; 28 } 29 } 30 int calc(int a,int b) 31 { 32 a=1ll*a*(a+1)/2%mod; 33 b=1ll*b*(b+1)/2%mod; 34 return (b-a+mod)%mod; 35 } 36 int find(int n) 37 { 38 if(n<=S)return miu[n]; 39 if(~a[nn/n])return a[nn/n]; 40 int ans=1,pos,sum; 41 for(int i=2;i<=n;i=pos+1) 42 { 43 pos=n/(n/i); 44 sum=1ll*calc(i-1,pos)*find(n/i)%mod; 45 ans=(ans-sum+mod)%mod; 46 } 47 return a[nn/n]=ans; 48 } 49 int solve(int a,int b) 50 { 51 a=find(a-1);b=find(b); 52 return (b-a+mod)%mod; 53 } 54 int work(int n) 55 { 56 int ans=0,pos; 57 for(int i=1;i<=n;i=pos+1) 58 { 59 pos=n/(n/i); 60 ans=(ans+1ll*calc(i-1,pos)*(n/i)%mod)%mod; 61 } 62 return 1ll*ans*ans%mod; 63 } 64 int main() 65 { 66 scanf("%d",&n); 67 pre();nn=n;memset(a,-1,sizeof(a)); 68 for(int i=1;i<=n;i=pos+1) 69 { 70 pos=n/(n/i); 71 ans=(ans+1ll*solve(i,pos)*work(n/i)%mod)%mod; 72 } 73 printf("%d",ans); 74 return 0; 75 }
9.【bzoj3512】DZY Loves Math IV
题意：见原题。
分析：permuiの博客
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<cmath> 5 #include<map> 6 #define LL long long 7 using namespace std; 8 const int N=1e6; 9 const int mod=1e9+7; 10 const int inv=(mod+1)/2; 11 int n,m,ans,tot,now; 12 int pri[N+5],phi[N+5],mn[N+5],a[N+5]; 13 bool np[N+5]; 14 map<int,int>M[N+5]; 15 int power(int a,int b) 16 { 17 int ans=1; 18 while(b) 19 { 20 if(b&1)ans=1ll*ans*a%mod; 21 a=1ll*a*a%mod;b>>=1; 22 } 23 return ans; 24 } 25 int getsum(int n) 26 { 27 if(n<=N)return phi[n]; 28 if(~a[m/n])return a[m/n]; 29 int ans=1ll*n*(n+1)%mod*inv%mod,pos; 30 for(int i=2;i<=n;i=pos+1) 31 { 32 pos=n/(n/i); 33 ans=(ans-1ll*(pos-i+1)*getsum(n/i)%mod+mod)%mod; 34 } 35 return a[m/n]=ans; 36 } 37 int Phi(int n) 38 { 39 int ans=1,num; 40 if(n<=N){ans=(phi[n]-phi[n-1]+mod)%mod;return ans;} 41 int m=sqrt(n)+0.5; 42 for(int i=1;pri[i]<=m;i++) 43 { 44 if(n%pri[i])continue; 45 num=0;while(n%pri[i]==0)num++,n/=pri[i]; 46 ans*=power(pri[i],num-1)*(pri[i]-1); 47 } 48 if(n!=1)ans*=(n-1); 49 return ans; 50 } 51 int S(int n,int m) 52 { 53 if(m==0)return 0; 54 if(n==1)return getsum(m); 55 if(M[n][m])return M[n][m]; 56 int ans=0; 57 for(int i=1;i*i<=n;i++) 58 { 59 if(n%i)continue; 60 int j=n/i;ans=(ans+1ll*Phi(j)*S(i,m/i)%mod)%mod; 61 if(i!=j)ans=(ans+1ll*Phi(i)*S(j,m/j)%mod)%mod; 62 } 63 return M[n][m]=ans; 64 } 65 int main() 66 { 67 phi[1]=mn[1]=1; 68 for(int i=2;i<=N;i++) 69 { 70 if(!np[i])pri[++tot]=i,phi[i]=i-1,mn[i]=i; 71 for(int j=1;i*pri[j]<=N;j++) 72 { 73 now=i*pri[j];np[now]=true; 74 if(i%pri[j]==0) 75 { 76 phi[now]=phi[i]*pri[j]; 77 mn[now]=mn[i];break; 78 } 79 phi[now]=phi[i]*(pri[j]-1); 80 mn[now]=mn[i]*pri[j]; 81 } 82 phi[i]=(phi[i-1]+phi[i])%mod; 83 } 84 memset(a,-1,sizeof(a)); 85 scanf("%d%d",&n,&m); 86 for(int i=1;i<=n;i++) 87 ans=(ans+1ll*i/mn[i]*S(mn[i],m)%mod)%mod; 88 printf("%d",ans); 89 return 0; 90 }
【洲阁筛】
〖相关资料〗
〖相关题目〗
1.【spoj】DIVCNT3 - Counting Divisors (cube)
题意：见原题。
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<cstdlib> 5 #include<algorithm> 6 #define LL long long 7 using namespace std; 8 const int N=320000; 9 const int mx=316240; 10 const int p=1000003; 11 int tot,id,pri[N],num[N],ci[N],mn[N]; 12 int cnt,first[p],w[p],nex[p],D[p]; 13 LL T,n,ans1,f[N],f2[p],sum[N],to[p],d[p],g[p]; 14 LL read() 15 { 16 LL x=0,f=1;char c=getchar(); 17 while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1;c=getchar();} 18 while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘;c=getchar();} 19 return x*f; 20 } 21 void ins(int x,LL y) 22 { 23 int id=y%p;to[++cnt]=y;w[cnt]=x; 24 nex[cnt]=first[id];first[id]=cnt; 25 } 26 int idx(LL x) 27 { 28 for(int i=first[x%p];i;i=nex[i]) 29 if(to[i]==x)return w[i]; 30 return 0; 31 } 32 void work1() 33 { 34 cnt=id=0;int k; 35 for(LL i=1;i<=n;i=n/(n/i)+1)first[n/i%p]=0; 36 for(LL i=1;i<=n;i=n/(n/i)+1) 37 ins(++id,n/i),d[id]=g[id]=n/i,D[id]=0; 38 for(int i=1;i<=tot;i++) 39 for(int j=1;j<=id&&1ll*pri[i]*pri[i]<=d[j];j++) 40 k=idx(d[j]/pri[i]),g[j]-=g[k]-(i-1-D[k]),D[j]=i; 41 } 42 void work2() 43 { 44 for(int i=tot;i>=1;i--) 45 for(int j=1;j<=id&&1ll*pri[i]*pri[i]<=d[j];j++) 46 { 47 if(pri[i+1]>d[j])f2[j]=1; 48 else if(1ll*pri[i+1]*pri[i+1]>d[j]) 49 f2[j]=(num[min(1ll*mx,d[j])]-num[pri[i+1]-1])*4+1; 50 for(LL pi=pri[i],c=1;pi<=d[j];pi*=pri[i],c++) 51 { 52 LL k=d[j]/pi,k2; 53 if(pri[i+1]>k)k2=1; 54 else if(1ll*pri[i+1]*pri[i+1]>k) 55 k2=(num[min(1ll*mx,k)]-num[pri[i+1]-1])*4+1; 56 else k2=f2[idx(k)]; 57 f2[j]+=k2*(3*c+1); 58 } 59 } 60 } 61 LL solve(LL n) 62 { 63 if(n<=mx)return sum[n]; 64 ans1=0;work1();work2(); 65 for(int i=1,pos;i<=mx;i=pos+1) 66 { 67 int j=idx(n/i);LL k; 68 if(pri[tot+1]>n/i)k=0; 69 else k=g[j]-(tot-D[j]+1); 70 ans1+=(sum[pos=min(1ll*mx,n/(n/i))]-sum[i-1])*k*4; 71 } 72 return ans1+f2[1]; 73 } 74 int main() 75 { 76 f[1]=sum[1]=1; 77 for(int i=2;i<=mx;i++) 78 { 79 num[i]=num[i-1]+(!f[i]); 80 if(!f[i])pri[++tot]=i,f[i]=4,mn[i]=i,ci[i]=1; 81 for(int j=1,now;(now=i*pri[j])<=mx;j++) 82 { 83 if(i%pri[j]==0) 84 { 85 ci[now]=ci[i]+1; 86 f[now]=f[i/mn[i]]*(ci[now]*3+1); 87 mn[now]=mn[i]*pri[j]; 88 break; 89 } 90 f[now]=f[i]*4;mn[now]=pri[j];ci[now]=1; 91 } 92 sum[i]=(sum[i-1]+f[i]); 93 } 94 pri[tot+1]=316241; 95 T=read(); 96 while(T--) 97 { 98 n=read(); 99 printf("%lld\\n",solve(n)); 100 } 101 return 0; 102 }
2.【51nod 1184】第N个质数
题意：给出一个数N，求第N个质数。
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #include<cmath> 5 #define LL long long 6 using namespace std; 7 const int N=12500000; 8 int n,tot,f[N+5],pri[N+5]; 9 LL l,r,mid,ans,tmp; 10 bool np[N+5]; 11 LL g(LL n,int m) 12 { 13 if(!m)return n; 14 if(m==1)return n-n/2; 15 if(n<=N) 16 { 17 if(m>=f[n])return 1; 18 if(m>=f[(int)sqrt(n)])return f[n]-m+1; 19 } 20 return g(n,m-1)-g(n/pri[m],m-1); 21 } 22 LL calc(LL n) 23 { 24 if(n<=N)return f[n]; 25 tmp=ceil(sqrt(n)); 26 return f[tmp]+g(n,f[tmp])-1; 27 } 28 int main() 29 { 30 for(int i=2;i<=N;i++) 31 { 32 if(!np[i])pri[++tot]=i; 33 for(int j=1;i*pri[j]<=N;j++) 34 { 35 np[i*pri[j]]=true; 36 if(i%pri[j]==0)break; 37 } 38 f[i]=f[i-1]+(!np[i]); 39 } 40 scanf("%d",&n); 41 l=1;r=22801763489; 42 if(n<=500000000)r=11037271757; 43 else l=11037271758; 44 if(n<=750000000)r=16875026921; 45 else l=16875026922; 46 if(n<=850000000)r=19236701629; 47 else l=19236701630; 48 if(n<=900000000)r=20422213579; 49 else l=20422213580; 50 if(n<=940000000)r=22801763489; 51 else l=22801763490; 52 if(n<=970000000)r=22086742277; 53 else l=22086742278; 54 if(n<=990000000)r=22563323957; 55 else l=22563323958; 56 while(l<=r) 57 { 58 mid=(l+r)>>1; 59 if(calc(mid)>=n)ans=mid,r=mid-1; 60 else l=mid+1; 61 } 62 printf("%lld",ans); 63 return 0; 64 }
以上是关于算法总结积性函数相关的主要内容，如果未能解决你的问题，请参考以下文章