【BZOJ 3512】 DZY Loves Math IV（杜教筛+记忆化搜索）

传送门

考虑$n$较小
设$S(n,m)=sum_{i=1}^mphi(in)$
设$n=prod_{i=1}^kp_i^{a_i},x=prod_{i=1}^kp_i,y=prod_{i=1}^kp_i^{a_i-1}$
则$n=x*y,phi(n)=phi(x)*y$
$S(n,m)=ysum_{i=1}^mphi(ix)$
设$d=gcd(x,i)$
$=ysum_{i=1}^mphi(i)phi(frac x d)d$
考虑用$Id=phi*I$代换
$=ysum_{i}phi(i)phi(frac x d)sum_{k|d}phi(k)$
$=ysum_{i}phi(i)sum_{k|d}phi(frac xk)$
$=ysum_{k|x}phi(frac x k)sum_{i=1}^{frac m k}phi(ki)$
$=ysum_{k|x}phi(frac x k)S(k,frac m k)$

先用杜教筛求出$S(1,m)$
然后利用这个式子记忆化搜索

复杂度据说是$O(nsqrt m+m^{frac 2 3})$
太菜完全不会算复杂度。。。。

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define ll long long
#define fi first
#define se second
#define bg begin
cs int RLEN=1<<20|1;
inline char gc(){
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob)?EOF:*ib++;
}
inline int read(){
    char ch=gc();
    int res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
cs int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?(a-mod):a;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return (r>=mod)?(r%mod):r;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
cs int N=100005,M=1000006;
int pr[M],phi[M],tot,X[N],mnp[N];
bitset<M> vis;
vector<int>pri[N];
inline void init(int len=M-6){
    phi[1]=1;
    for(int i=2;i<=len;i++){
        if(!vis[i])pr[++tot]=i,phi[i]=i-1;
        for(int j=1,p;j<=tot&&i*pr[j]<=len;j++){
            p=i*pr[j],vis[p]=1;
            if(i%pr[j]==0){
                phi[p]=phi[i]*pr[j];break;
            }
            phi[p]=phi[i]*phi[pr[j]];
        }
    }
    for(int i=1;i<=len;i++)Add(phi[i],phi[i-1]);
    len=N-5;
    for(int i=1;i<=len;i++)pri[i].pb(1);
    for(int i=2;i<=len;i++)
    for(int j=i;j<=len;j+=i){
        if(!mnp[j])mnp[j]=i;
        pri[j].pb(i);
    }
    for(int i=1;i<=len;i++){
        int x=i;X[i]=1;
        while(x>1){if(X[i]%mnp[x])X[i]*=mnp[x];x/=mnp[x];}
    }
}
inline int S(int x){return 1ll*x*(x+1)/2%mod;}
map<int,int>mp;
inline int calc(int n){
    if(n<=M-6)return phi[n];
    if(mp.count(n))return mp[n];
    int ret=S(n);
    for(int i=2,j;i<=n;i=j+1){
        j=(n/(n/i));
        Dec(ret,mul(j-i+1,calc(n/i)));
    }
    return mp[n]=ret;
}
map<int,int> s[N];
inline int solve(int n,int m){
    if(n==1)return calc(m);
    if(m==1)return dec(phi[n],phi[n-1]);
    if(!m)return 0;
    if(s[n].count(m))return s[n][m];
    int x=X[n],y=n/x;
    int ret=0;
    for(int i=0,k;i<pri[x].size();i++){
        k=pri[x][i];
        Add(ret,mul(dec(phi[x/k],phi[x/k-1]),solve(k,m/k)));
    }
    Mul(ret,y);
    s[n][m]=ret;return ret;
}
int n,m;
int main(){
    #ifdef Stargazer
    freopen("lx.in","r",stdin);
    #endif
    init();
    n=read(),m=read();
    calc(m);
    int ret=0;
    for(int i=1;i<=n;i++)Add(ret,solve(i,m));
    cout<<ret;
}