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

题目链接
题意：求n<=100000,m<=1000000000。
题解：首先我们有一个结论。
若)。
怎么证？
)
推导过程中用到了很多|=1的性质。
这个式子就很显然了，因为根据)，想一想就知道了。于是就证完了！
若)=1。
则)∗nk
接下来我们继续推导。
我们令)
则)。
对于)
否则若n≠1，
设k为最小的正整数满足)∗nk
否则当)
跟求)=i。
=>)2
=>)2
=>i=1带入。
=>)。
于是我们只需要开个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
",(ans+mod)%mod);
    return 0;
}