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

题目链接
题意：求∑ni=1∑mj=1φ(ij)$\sum _{i=1}^n\sum _{j=1}^m\phi (ij)$模1000000007$1000000007$的值。n<=100000,m<=1000000000$n<=100000,m<=1000000000$。
题解：首先我们有一个结论。
若|μ(n)|=1$|\mu (n)|=1$，则φ(nm)=∑k|gcd(n,m)φ(nk)φ(m)$\phi (nm)=\sum _{k|gcd(n,m)}\phi (\frac{n}{k})\phi (m)$。
怎么证？
∑k|gcd(n,m)φ(nk)φ(m)$\sum _{k|gcd(n,m)}\phi (\frac{n}{k})\phi (m)$
=φ(m)∑k|gcd(n,m)φ(nk)$=\phi (m)\sum _{k|gcd(n,m)}\phi (\frac{n}{k})$
=φ(m)∑k|gcd(n,m)φ(n)φ(k)$=\phi (m)\sum _{k|gcd(n,m)}\frac{\phi (n)}{\phi (k)}$
=φ(n)φ(m)∑k|gcd(n,m)1φ(k)$=\phi (n)\phi (m)\sum _{k|gcd(n,m)}\frac{1}{\phi (k)}$
=φ(n)φ(m)∑k|gcd(n,m)φ(gcd(n,m))φ(k)φ(gcd(n,m))$=\phi (n)\phi (m)\sum _{k|gcd(n,m)}\frac{\frac{\phi (gcd(n,m))}{\phi (k)}}{\phi (gcd(n,m))}$
=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m))φ(k)$=\phi (n)\phi (m)\frac{1}{\phi (gcd(n,m))}\sum _{k|gcd(n,m)}\frac{\phi (gcd(n,m))}{\phi (k)}$
=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m)k)$=\phi (n)\phi (m)\frac{1}{\phi (gcd(n,m))}\sum _{k|gcd(n,m)}\phi (\frac{gcd(n,m)}{k})$
=φ(n)φ(m)1φ(gcd(n,m))∑k|gcd(n,m)φ(gcd(n,m)k)$=\phi (n)\phi (m)\frac{1}{\phi (gcd(n,m))}\sum _{k|gcd(n,m)}\phi (\frac{gcd(n,m)}{k})$
=φ(n)φ(m)gcd(n,m)φ(gcd(n,m))$=\phi (n)\phi (m)\frac{gcd(n,m)}{\phi (gcd(n,m))}$ 这是因为∑d|nφ(d)=n$\sum _{d|n}\phi (d)=n$。
=φ(ngcd(n,m))φ(m)gcd(n,m)$=\phi (\frac{n}{gcd(n,m)})\phi (m)gcd(n,m)$
推导过程中用到了很多|μ(n)|=1$|\mu (n)|=1$的性质。
这个式子就很显然了，因为根据|μ(n)|=1$|\mu (n)|=1$，ngcd(n,m)$\frac{n}{gcd(n,m)}$与m$m$互质，于是φ(ngcd(n,m))φ(m)gcd(n,m)=φ(nm)$\phi (\frac{n}{gcd(n,m)})\phi (m)gcd(n,m)=\phi (nm)$，想一想就知道了。于是就证完了！
若|μ(n)|=0$|\mu (n)|=0$，设k为最小的正整数满足k|n$k|n$，μ(k)=1$\mu (k)=1$。
则φ(nm)=φ(km)∗nk$\phi (nm)=\phi (km)\ast \frac{n}{k}$
接下来我们继续推导。
我们令S(n,m)=∑mi=1φ(ni)$S(n,m)=\sum _{i=1}^m\phi (ni)$
则ans=∑ni=1S(i,m)$ans=\sum _{i=1}^{n}S(i,m)$。
对于S(n,m)$S(n,m)$，若|μ(n)|=1$|\mu (n)|=1$且n≠1$n\ne 1$，则
S(n,m)$S(n,m)$
=∑mi=1φ(ni)$=\sum _{i=1}^m\phi (ni)$
=∑mi=1∑d|gcd(n,i)φ(nd)φ(i)$=\sum _{i=1}^m\sum _{d|gcd(n,i)}\phi (\frac{n}{d})\phi (i)$
=∑d|nφ(nd)∑⌊md⌋i=1φ(i)$=\sum _{d|n}\phi (\frac{n}{d})\sum _{i=1}^{\lfloor \frac{m}{d}\rfloor }\phi (i)$
=∑d|nφ(nd)S(d,⌊md⌋)$=\sum _{d|n}\phi (\frac{n}{d})S(d,\lfloor \frac{m}{d}\rfloor )$
否则若|μ(n)|=0$|\mu (n)|=0$且n≠1$n\ne 1$，
设k为最小的正整数满足k|n$k|n$，μ(k)=1$\mu (k)=1$，则
S(n,m)=S(k,m)∗nk$S(n,m)=S(k,m)\ast \frac{n}{k}$
否则当n=1$n=1$，
S(n,m)=∑mi=1φ(i)$S(n,m)=\sum _{i=1}^m\phi (i)$
跟求μ$\mu$的前缀和类似，∑ni=1∑j|iφ(j)=n(n+1)2$\sum _{i=1}^n\sum _{j|i}\phi (j)=\frac{n(n+1)}{2}$，因为∑j|iφ(j)=i$\sum _{j|i}\phi (j)=i$。
=>∑nj=1∑⌊nj⌋i=1φ(j)=n(n+1)2$\sum _{j=1}^n\sum _{i=1}^{\lfloor \frac{n}{j}\rfloor }\phi (j)=\frac{n(n+1)}{2}$
=>∑ni=1∑⌊ni⌋j=1φ(j)=n(n+1)2$\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\phi (j)=\frac{n(n+1)}{2}$
=>∑nj=1φ(j)=n(n+1)2−∑ni=2∑⌊ni⌋j=1φ(j)$\sum _{j=1}^n\phi (j)=\frac{n(n+1)}{2}-\sum _{i=2}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }\phi (j)$，即把i=1$i=1$带入。
=>S(n,m)=m(m+1)2−∑mi=2S(n,⌊mi⌋)$S(n,m)=\frac{m(m+1)}{2}-\sum _{i=2}^{m}S(n,\lfloor \frac{m}{i}\rfloor )$。
于是我们只需要开个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;
}