【bzoj3512】DZY Loves Math IV 【杜教筛】

题目链接
题意：求∑ni=1∑mj=1φ(ij)∑i=1n∑j=1mφ(ij)模10000000071000000007的值。n<=100000,m<=1000000000n<=100000,m<=1000000000。
题解：首先我们有一个结论。
若|μ(n)|=1|μ(n)|=1，则φ(nm)=∑k|gcd(n,m)φ(nk)φ(m)φ(nm)=∑k|gcd(n,m)φ(nk)φ(m)。
怎么证？
∑k|gcd(n,m)φ(nk)φ(m)∑k|gcd(n,m)φ(nk)φ(m)
=φ(m)∑k|gcd(n,m)φ(nk)=φ(m)∑k|gcd(n,m)φ(nk)
=φ(m)∑k|gcd(n,m)φ(n)φ(k)=φ(m)∑k|gcd(n,m)φ(n)φ(k)
=φ(n)φ(m)∑k|gcd(n,m)1φ(k)=φ(n)φ(m)∑k|gcd(n,m)1φ(k)
=φ(n)φ(m)∑k|gcd(n,m)φ(gcd(n,m))φ(k)φ(gcd(n,m))=φ(n)φ(m)∑k|gcd(n,m)φ(gcd(n,m))φ(k)φ(gcd(n,m))
=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m))φ(k)=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m))φ(k)
=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m)k)=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m)k)
=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m)k)=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m)k)
=φ(n)φ(m)gcd(n,m)φ(gcd(n,m))=φ(n)φ(m)gcd(n,m)φ(gcd(n,m)) 这是因为∑d|nφ(d)=n∑d|nφ(d)=n。
=φ(ngcd(n,m))φ(m)gcd(n,m)=φ(ngcd(n,m))φ(m)gcd(n,m)
推导过程中用到了很多|μ(n)|=1|μ(n)|=1的性质。
这个式子就很显然了，因为根据|μ(n)|=1|μ(n)|=1，ngcd(n,m)ngcd(n,m)与mm互质，于是φ(ngcd(n,m))φ(m)gcd(n,m)=φ(nm)φ(ngcd(n,m))φ(m)gcd(n,m)=φ(nm)，想一想就知道了。于是就证完了！
若|μ(n)|=0|μ(n)|=0，设k为最小的正整数满足k|nk|n，μ(k)=1μ(k)=1。
则φ(nm)=φ(km)∗nkφ(nm)=φ(km)∗nk
接下来我们继续推导。
我们令S(n,m)=∑mi=1φ(ni)S(n,m)=∑i=1mφ(ni)
则ans=∑ni=1S(i,m)ans=∑i=1nS(i,m)。
对于S(n,m)S(n,m)，若|μ(n)|=1|μ(n)|=1且n≠1n≠1，则
S(n,m)S(n,m)
=∑mi=1φ(ni)=∑i=1mφ(ni)
=∑mi=1∑d|gcd(n,i)φ(nd)φ(i)=∑i=1m∑d|gcd(n,i)φ(nd)φ(i)
=∑d|nφ(nd)∑⌊md⌋i=1φ(i)=∑d|nφ(nd)∑i=1⌊md⌋φ(i)
=∑d|nφ(nd)S(d,⌊md⌋)=∑d|nφ(nd)S(d,⌊md⌋)
否则若|μ(n)|=0|μ(n)|=0且n≠1n≠1，
设k为最小的正整数满足k|nk|n，μ(k)=1μ(k)=1，则
S(n,m)=S(k,m)∗nkS(n,m)=S(k,m)∗nk
否则当n=1n=1，
S(n,m)=∑mi=1φ(i)S(n,m)=∑i=1mφ(i)
跟求μμ的前缀和类似，∑ni=1∑j|iφ(j)=n(n+1)2∑i=1n∑j|iφ(j)=n(n+1)2，因为∑j|iφ(j)=i∑j|iφ(j)=i。
=>∑nj=1∑⌊nj⌋i=1φ(j)=n(n+1)2∑j=1n∑i=1⌊nj⌋φ(j)=n(n+1)2
=>∑ni=1∑⌊ni⌋j=1φ(j)=n(n+1)2∑i=1n∑j=1⌊ni⌋φ(j)=n(n+1)2
=>∑nj=1φ(j)=n(n+1)2−∑ni=2∑⌊ni⌋j=1φ(j)∑j=1nφ(j)=n(n+1)2−∑i=2n∑j=1⌊ni⌋φ(j)，即把i=1i=1带入。
=>S(n,m)=m(m+1)2−∑mi=2S(n,⌊mi⌋)S(n,m)=m(m+1)2−∑i=2mS(n,⌊mi⌋)。
于是我们只需要开个map无脑记搜乱搞即可。时间复杂度不会算= =
丑得不堪入目的 代码：

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<map>
using namespace std;
typedef long long ll;
const int N=1000005;
const ll mod=1000000007;
int n,m,p[N];
bool vis[N];
ll ans,mu[N],phi[N],sum[N];
map<int,map<int,ll> >mp;
ll solve(int n,int m){
    if(m<=1){
        return phi[n*m];
    }else if(n==1){
        if(m<=1000000){
            return sum[m];
        }
        if(mp[n][m]){
            return mp[n][m];
        }
        ll res=1LL*m*(m+1)/2;
        for(int i=2,last;i<=m;i=last+1){
            last=m/(m/i);
            res-=solve(n,m/i)*(last-i+1)%mod;
            res%=mod;
        }
        return mp[n][m]=res;
    }else{
        if(mp[n][m]){
            return mp[n][m];
        }
        int tmp=0;
        for(int i=1;i*i<=n;i++){
            if(n%i==0&&mu[n/i]){
                tmp=i;
            break;
            }
        }
        if(!tmp){
            for(int i=sqrt(n);i>=1;i--){
                if(n%i==0&&mu[i]){
                    tmp=n/i;
                    break;
                }
            }
        }
        n/=tmp;
        ll res=0;
        for(int i=1;i*i<=n;i++){
            if(n%i==0){
                res+=phi[n/i]*solve(i,m/i)%mod;
                res%=mod;
                if(i*i!=n){
                    res+=phi[i]*(solve(n/i,m/(n/i)))%mod;
                    res%=mod;
                }
            }
        }
        n*=tmp;
        res=res*tmp%mod;
        return mp[n][m]=res;
    }
}
int main(){
    mu[1]=phi[1]=1;
    for(int i=2;i<=1000000;i++){
        if(!vis[i]){
            p[++p[0]]=i;
            mu[i]=-1;
            phi[i]=i-1;
        }
        for(int j=1;j<=p[0]&&i*p[j]<=1000000;j++){
            vis[i*p[j]]=true;
            if(i%p[j]){
                mu[i*p[j]]=-mu[i];
                phi[i*p[j]]=phi[i]*(p[j]-1);
            }else{
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
        }
    }
    for(int i=1;i<=1000000;i++){
        sum[i]=sum[i-1]+phi[i];
        sum[i]%=mod;
    }
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;i++){
        ans+=solve(i,m);
        ans%=mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}

转载于:https://www.cnblogs.com/2016gdgzoi471/p/9476863.html