1. 程式人生 > >bzoj 4916: 神犇和蒟蒻【歐拉函數+莫比烏斯函數+杜教篩】

bzoj 4916: 神犇和蒟蒻【歐拉函數+莫比烏斯函數+杜教篩】

pac pan using 函數 莫比烏斯函數 right for body ace

居然扒到了學長出的題
和3944差不多(?),雖然一眼看上去很可怕但是仔細觀察發現,對於mu來講,答案永遠是1(對於帶平方的,mu值為0,1除外),然後根據歐拉篩的原理,\( \sum_{i=1}^{n}\phi(i^2)=\sum_{i=1}^{n}\phi(i)*i \),然後就可以正常推了:

\[ g(n)=\sum_{i=1}^{n}i\sum_{d=1}^{i}[d|i]\phi(d)=\sum_{i=1}^{n}i^2=\frac{n(n+1){2n+1}}{6} \]
\[ s(n)=\sum_{i=1}^{n}i\phi(i) \]那麽把g展開:
\[ g(n)=\sum_{i=2}^{n}i\sum_{d=1}^{i-1}[d|i]\phi(d)+s(n) \]


\[ s(n)=g(n)-\sum_{i=2}^{n}i\sum_{d=1}^{i-1}[d|i]\phi(d) \]
\[ =g(n)-\sum_{k=2}^{n}k\sum_{d=1}^{\left \lfloor \frac{n}{k} \right \rfloor}d\phi(d) \]
\[ =g(n)-\sum_{k=2}^{n}k*s(\left \lfloor \frac{n}{k} \right \rfloor) \]
\[ =\frac{n(n+1){2n+1}}{6}-\sum_{k=2}^{n}k*s(\left \lfloor \frac{n}{k} \right \rfloor) \]

這就是標準的杜教篩遞歸形式了,直接求解即可。

#include<iostream>
#include<cstdio>
#include<map>
#include<cstring>
using namespace std;
const int N=1000005,m=1000000,mod=1e9+7;
int phi[N],mi[N],q[N],tot,n,k,s[N],ans[N];
bool v[N];
map<long long,int>mp;
int S(int n,int l)
{
    if(l<=1)
        return
phi[n*l]; if(n==1) { if(l<=m) return s[l]; if(ans[k/l]!=-1) return ans[k/l]; long long re=(long long)l*(l+1)/2%mod; for(int i=2,la;i<=l;i=la+1) { la=l/(l/i); if(l/i<=m) re=(re-(long long)s[l/i]*(la-i+1)%mod)%mod; else re=(re-(long long)S(1,l/i)*(la-i+1)%mod)%mod; } return ans[k/l]=(re%mod+mod)%mod; } if(mp[(long long)n*mod+l]) return mp[(long long)n*mod+l]; long long re=0ll; for(int i=1;i*i<=n;i++) if(n%i==0) { re=(re+(long long)phi[n/i]*S(i,l/i)%mod)%mod; if(i*i!=n) re=re+(long long)phi[i]*S(n/i,l/(n/i))%mod; } return mp[(long long)n*mod+l]=(re%mod+mod)%mod; } // int S(int n,int l) // { // if (l<=1) return phi[n*l]; // if (n==1) // { // if (l<=m) return s[l]; // if (ans[k/l]!=-1) return ans[k/l]; // int re=(int)l*(l+1)/2%mod; // for (int i=2,la;i<=l;i=la+1) // { // la=l/(l/i); // if (l/i<=m) re=re-(int)(la-i+1)*s[l/i]%mod+mod; // else re=re-(int)(la-i+1)*S(1,l/i)%mod+mod; // } // return ans[k/l]=re%mod; // } // else // { // if (mp[(int)n*mod+l]) return mp[(int)n*mod+l]; // int re=0; // for (int i=1;i*i<=n;i++) // if (n%i==0) // { // re=re+(int)phi[n/i]*S(i,l/i)%mod; // if (i*i!=n) re=re+(int)phi[i]*S(n/i,l/(n/i))%mod; // } // return mp[(int)n*mod+l]=re%mod; // } // } int main() { memset(ans,-1,sizeof(ans)); mi[1]=phi[1]=1; for(int i=2;i<=m;i++) { if(!v[i]) { q[++tot]=i; phi[i]=i-1; mi[i]=i; } for(int j=1;j<=tot&&i*q[j]<=m;j++) { int k=i*q[j]; v[k]=1; if(i%q[j]==0) { phi[k]=phi[i]*q[j]; mi[k]=mi[i]; break; } phi[k]=phi[i]*(q[j]-1); mi[k]=mi[i]*q[j]; } } for(int i=1;i<=m;i++) s[i]=(s[i-1]+phi[i])%mod; scanf("%lld%lld",&n,&k); if(n>k) swap(n,k); long long ans=0ll; for(int i=1;i<=n;i++) ans=(ans+((long long)i/mi[i]*S(mi[i],k)%mod))%mod; printf("%lld\n",(ans%mod+mod)%mod); return 0; }

bzoj 4916: 神犇和蒟蒻【歐拉函數+莫比烏斯函數+杜教篩】