bzoj3512 DZY Loves Math IV

因为n比较小，m比较大，我们枚举n，计算

s(n,m)=∑i=1mφ(ni)

首先把n的所有完全平方因子提出来，也就是使得|μ(n)|=1。接下来任取n的一个质因子p，可以发现当i|p时φ(ni)=pφ(nip)，否则φ(ni)=(p−1)φ(nip)。因此

s(n,m)=(p−1)s(np,m)+s(n,⌊mp⌋)

当m=1时，答案就是φ(n)。当n=1时，用杜教筛求解。

#include<cstdio>#include<cmath>#include<algorithm>using namespace std;#define LL long longconst int maxn=1000000,p=1000000007;int dv[maxn+10],phi[maxn+10],sphi[maxn+10],prm[maxn+10],vis[maxn+10],tot;inline int inc(int x,int y){    x+=y;    return x>=p?x-p:x;}inline int dec(int x,int y){    x-=y;    return x<0?x+p:x;}struct one_hash{    int fir[maxn+10],ne[maxn+10],val[maxn+10],ans[maxn+10],tot,p;    void init()    {        p=1000007;    }    inline int find(int n)    {        for (int i=fir[n%p];i;i=ne[i])            if (val[i]==n) return ans[i];        return -1;    }    inline void ins(int n,int x)    {        if (tot>maxn+5) return;        int y=n%p;        tot++;        ne[tot]=fir[y];        fir[y]=tot;        val[tot]=n;        ans[tot]=x;    }}h1;void init(){    phi[1]=sphi[1]=1;    for (int i=2;i<=maxn;i++)    {        if (!dv[i])        {            phi[i]=i-1;            prm[++tot]=i;            dv[i]=i;        }        for (int j=1;j<=tot&&(LL)i*prm[j]<=maxn;j++)            if (i%prm[j])            {                dv[i*prm[j]]=prm[j];                phi[i*prm[j]]=phi[i]*(prm[j]-1);            }            else            {                dv[i*prm[j]]=dv[i]*prm[j];                phi[i*prm[j]]=phi[i]*prm[j];                break;            }        sphi[i]=inc(sphi[i-1],phi[i]);    }    h1.init();}inline int sum(int n){    if (n<=maxn) return sphi[n];    int ret=h1.find(n);    if (ret!=-1) return ret;    ret=(LL)n*(n+1)/2%p;    for (int i=2,j;i<=n;i=j+1)    {        j=n/(n/i);        ret=dec(ret,(LL)(j-i+1)*sum(n/i)%p);    }    h1.ins(n,ret);    return ret;}int get(int n){    int m=sqrt(n+0.5),ret=1;    for (int i=1;i<=tot&&prm[i]<=m;i++)        if (n%prm[i]==0)        {            ret*=prm[i];            while (n%prm[i]==0) n/=prm[i];        }    if (n>1) ret*=n;    return ret;}inline int solve(int n,int m){    if (n==1) return sum(m);    if (m==1) return phi[n];    if (n==0||m==0) return 0;    return inc((LL)(dv[n]-1)*solve(n/dv[n],m)%p,solve(n,m/dv[n]));}int main(){    int n,m,x,ans=0;    init();    scanf("%d%d",&n,&m);    for (int i=1;i<=n;i++) vis[i]=-1;    for (int i=1;i<=n;i++)    {        x=get(i);        if (vis[x]==-1) vis[x]=solve(x,m);        ans=inc(ans,(LL)i/x*vis[x]%p);    }    printf("%d\n",ans);}