BZOJ 3512: DZY Loves Math IV [杜教筛]

BZOJ 3512: DZY Loves Math IV [杜教筛]

注意是去除所有质因子的乘积的n

厉害的思想是构造互质情况把phi(i*j)分开

补充candy？的最后一步推导：枚举e，那么gcd(i,n)要是e的倍数，枚举是e的多少倍（上界[m/e]），乘上phi(i*d)，就是S(n,m)的形式了

#include<bits/stdc++.h>
#define reg register int
#define il inline
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
namespace Miracle{
const int N=1e5+5;
const int M=1999880;
const int mod=1e9+7;
int pri[M],tot;
int pro[M];
int phi[M],vis[M];
int sum[M];
int n,m;
int ad(int x,int y){
    return (x+y)>=mod?x+y-mod:x+y;
}
void sieve(int n){
    phi[1]=sum[1]=1;
    pro[1]=1;
    for(reg i=2;i<=n;++i){
        if(!vis[i]){
            pri[++tot]=i;
            phi[i]=i-1;
            pro[i]=i;
        }
        for(reg j=1;j<=tot;++j){
            if(i*pri[j]>n) break;
            vis[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                pro[i*pri[j]]=pro[i];
                break;
            }
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
            pro[i*pri[j]]=pro[i]*pri[j];
        }
        sum[i]=ad(sum[i-1],phi[i]);
    }
}
const int G=31640;
int p1[G],p2[G];
int Phi(int n){
    if(n<=M-5){
        return sum[n];
    }
    //cout<<" 23333 ------"<<endl;
    if(n<G){
        if(p1[n]!=-1) return p1[n];
    }
    else{
        if(p2[m/n]!=-1) return p2[m/n];
    }
    int ret=(ll)n*(n+1)/2%mod;
    for(reg i=2,x=0;i<=n;i=x+1){
        x=(n/(n/i));
        ret=ad(ret,mod-1LL*(x-i+1)*Phi(n/i)%mod);
    }
    if(n<G){
        return p1[n]=ret;
    }else{
        return p2[m/n]=ret;
    }
}
map<int,int>mp[N];
int S(int n,int m){
    //cout<<" S "<<n<<" "<<m<<endl;
    if(m==0||n==0) return 0;
    if(n==1) return Phi(m);
    if(m==1) return phi[n];
    if(mp[n][m]) return mp[n][m];
    int ret=0;
    for(reg i=1;i*i<=n;++i){
        if(n%i==0){
            ret=ad(ret,(ll)phi[n/i]*S(i,m/i)%mod);
            if(i!=n/i) ret=ad(ret,(ll)phi[i]*S(n/i,m/(n/i))%mod);
        }
    }
    return mp[n][m]=ret;
}
int main(){
    rd(n);rd(m);
    if(n>m) swap(n,m);
    sieve(M-5);
//    cout<<" after sieve "<<endl;
    ll ans=0;
    memset(p1,-1,sizeof p1);
    memset(p2,-1,sizeof p2);
    for(reg i=1;i<=n;++i) {
    //    cout<<" ii "<<i<<endl;
        ans=(ans+(ll)i/pro[i]*S(pro[i],m)%mod)%mod;
    }
    printf("%lld",ans);
    return 0;
}

}
signed main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2019/3/7 20:51:51
*/