bzoj4916 神犇和蒟蒻

来源:互联网 发布:网络平台合作协议 编辑:程序博客网 时间:2024/06/04 17:49

μ的不用说了。
构造和式

i=1nd|iφ(d)did

一方面=ni=1i2
另一方面,=nx=1xnxd=1φ(d)d
而我们又知道φ(i2)=iφ(i)
因此不难杜教筛了。

#include<cstdio>#include<algorithm>using namespace std;#define LL long longconst int p=1000000007,mod=1000007,maxn=1000000,inv=166666668;int vis[maxn],prm[maxn],phi[maxn],sum[maxn];struct hash{    int fir[maxn+10],ne[maxn+10],val[maxn+10],ans[maxn+10],tot;    int find(int n)    {        for (int i=fir[n%mod];i;i=ne[i])            if (val[i]==n) return ans[i];        return -1;    }    void ins(int n,int x)    {        int y=n%mod;        tot++;        ne[tot]=fir[y];        fir[y]=tot;        val[tot]=n;        ans[tot]=x;    }}h;int inc(int x,int y){    x+=y;    return x>p?x-p:x;}int dec(int x,int y){    x-=y;    return x<0?x+p:x;}int solve(int n){    if (n<=maxn) return sum[n];    int ret=h.find(n);    if (ret!=-1) return ret;    ret=(LL)n*(n+1)%p*(2*n+1)%p*inv%p;    for (int i=2,j;i<=n;i=j+1)    {        j=n/(n/i);        ret=dec(ret,(LL)(i+j)*(j-i+1)/2*solve(n/i)%p);    }    h.ins(n,ret);    return ret;}int main(){    int tot=0,n;    sum[1]=1;    for (int i=2;i<=maxn;i++)    {        if (!vis[i])        {            prm[++tot]=i;            phi[i]=i-1;        }        for (int j=1;j<=tot&&(LL)i*prm[j]<=maxn;j++)        {            vis[i*prm[j]]=1;            if (i%prm[j]) phi[i*prm[j]]=phi[i]*(prm[j]-1);            else            {                phi[i*prm[j]]=phi[i]*prm[j];                break;            }        }        sum[i]=inc(sum[i-1],(LL)i*phi[i]%p);    }    scanf("%d",&n);    printf("1\n%d\n",solve(n));}
原创粉丝点击