【算法总结】积性函数相关

【线性筛】

〖模板代码〗

  [线性筛质数]

int main()
{ 
    n=read();
    for(int i=2;i<=n;i++)
    {
        if(!f[i])pri[++cnt]=i;
        for(int j=1;j<=cnt;j++)
        {
            if(i*pri[j]>n)break;
            f[i*pri[j]]=1;
            if(i%pri[j]==0)break;
        }
    }
}

  [线性求约数和]

void pre()
{
    miu[1]=1;f[1].first=1;
    for(int i=2;i<=mx;i++)
    {
        if(!np[i]){pri[++cnt]=i;d[i]=i;f[i].first=i+1;miu[i]=-1;}
        for(int j=1;i*pri[j]<=mx;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0)
            {
                miu[i*pri[j]]=0;
                d[i*pri[j]]=d[i]*pri[j];
                f[i*pri[j]].first=f[i].first+f[i/d[i]].first*d[i*pri[j]];
                break;
            }
            miu[i*pri[j]]=-miu[i];
            d[i*pri[j]]=pri[j];
            f[i*pri[j]].first=f[i].first*f[pri[j]].first;
        }
    }
    for(int i=1;i<=mx;i++)f[i].second=i;
}

〖相关题目〗

1.【Technocup 2018 - Elimination Round 2】F. Paths

题意：给定数字n，建立一个无向图。对于所有1~n之间的数字，当数字gcd(u,v)≠1时将u、v连一条边，边权为1。d(u, v)表示u到v的最短路，求所有d(u, v)的和，其中1 ≤ u < v ≤ n。

分析：Zsnuoの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
using namespace std;
const int N=1e7+5;
int n,m,tot,now,pri[N],p[N],phi[N],num[N],sum[N]; 
LL one,two,three;
int main()
{
    scanf("%d",&n);
    phi[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!p[i]){p[i]=pri[++tot]=i;phi[i]=i-1;}
        for(int j=1;j<=tot;j++)
        {
            if(i*pri[j]>n)break;
            p[i*pri[j]]=pri[j];
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
            else phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
    } 
    for(int i=2;i<=n;i++)one+=i-1-phi[i];
    for(int i=2;i<=n;i++)num[p[i]]++;
    for(int i=2;i<=n;i++)sum[i]=sum[i-1]+num[i];
    for(int i=2;i<=n;i++)two+=1ll*num[i]*sum[n/i];
    for(int i=2;i<=n;i++)if(1ll*p[i]*p[i]<=n)two--;
    two=two/2-one;m=n-1;
    for(int i=tot;i>=1;i--)
        if(pri[i]*2>n)m--;
        else break;
    three=1ll*m*(m-1)/2-one-two;
    printf("%I64d\n",one+two*2+three*3);
    return 0;
}

2.【51nod1584】加权约数和

题意：见原题

分析：KsClaの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=1e6;
const int mod=1e9+7;
int T,n,tot,now,pri[N+5],miu[N+5];
int mn[N+5],h[N+5],p[N+5],g[N+5],f[N+5];
bool np[N+5];
int read()
{
    int x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int main()
{
    h[1]=p[1]=miu[1]=g[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i])
        {
            pri[++tot]=i;
            miu[i]=-1;mn[i]=i;h[i]=i+1;
            p[i]=(1ll*i*i%mod+i+1)%mod;
        }
        for(int j=1;i*pri[j]<=N;j++)
        {
            now=i*pri[j];np[now]=true;
            if(i%pri[j]==0)
            {
                miu[now]=0;mn[now]=mn[i]*pri[j];
                h[now]=(h[i]+1ll*h[i/mn[i]]*mn[now]%mod)%mod;
                p[now]=(1ll*(mn[now]+mn[i])%mod*mn[now]%mod*p[i/mn[i]]%mod+p[i])%mod;
                break;
            }
            miu[now]=-miu[i];mn[now]=pri[j];
            h[now]=1ll*h[i]*h[pri[j]]%mod;
            p[now]=1ll*p[i]*p[pri[j]]%mod;
        }
    }
    for(int i=2;i<=N;i++)p[i]=1ll*p[i]*i%mod;
    for(int i=2;i<=N;i++)p[i]=(p[i]+p[i-1])%mod;
    for(int i=2;i<=N;i++)miu[i]=1ll*(miu[i]+mod)%mod*i%mod*i%mod;
    for(int i=2;i<=N;i++)g[i]=(g[i-1]+h[i])%mod;
    for(int i=2;i<=N;i++)g[i]=1ll*i*g[i]%mod*h[i]%mod;
    for(int i=1;i<=N;i++)
        for(int j=1;i*j<=N;j++)
            f[i*j]=(f[i*j]+1ll*miu[i]*g[j]%mod)%mod;
    for(int i=2;i<=N;i++)f[i]=(f[i-1]+f[i])%mod;
    T=read();
    for(int i=1;i<=T;i++)
    {
        n=read();
        n=(2*f[n]%mod-p[n]+mod)%mod;
        printf("Case #%d: %d\n",i,n);
    }
    return 0;
}

【杜教筛】

〖相关资料〗

《浅谈一类积性函数的前缀和》        《论逗逼的自我修养之寒假颓废记》

〖相关题目〗

1.【51nod1239】欧拉函数之和

题意：求φ的前缀和。

分析：ONION_CYCの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=5e6;
const int mod=1e9+7;
const int inv=(mod+1)/2;
int s,tot,pri[N+5],phi[N+5],a[N+5],b[N+5];
bool isp[N+5];
LL n;
int memory(LL x)
{
    if(x<=s)return a[x];
    else return b[n/x];
}
void add(LL x,int ans)
{
    if(x<=s)a[x]=ans;
    else b[n/x]=ans;
}
int solve(LL n)
{
    if(n<=N)return phi[n];
    if(~(memory(n)))return memory(n);
    int ans=1ll*(n%mod)*((n+1)%mod)%mod*inv%mod;
    LL pos;
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans-1ll*(pos-i+1)%mod*solve(n/i)%mod+mod)%mod;
    }
    add(n,ans);
    return ans;
}
int main()
{
    memset(a,-1,sizeof(a));
    memset(b,-1,sizeof(b));
    scanf("%lld",&n);
    s=sqrt(n);phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!isp[i]){pri[++tot]=i;phi[i]=i-1;}
        for(int j=1;j<=tot;j++)
        {
            if(i*pri[j]>N)break;
            isp[i*pri[j]]=true;
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
            else phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
        phi[i]=(phi[i]+phi[i-1])%mod;
    }
    if(n<=N)printf("%d",phi[n]);
    else printf("%d",solve(n));
    return 0;
}

2.【51nod1244】莫比乌斯函数之和

题意：求μ的前缀和。

分析：无

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=5e6;
int s,tot,pri[N+5],miu[N+5],a[N+5];
bool np[N+5];
LL nn,L,R;
LL solve(LL n)
{
    if(n<=N)return miu[n];
    if(~a[nn/n])return a[nn/n];
    LL ans=1,pos;
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=ans-(pos-i+1)*solve(n/i);
    }
    return a[nn/n]=ans;
}
int main()
{
    miu[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i]){miu[i]=-1;pri[++tot]=i;}
        for(int j=1;i*pri[j]<=N;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0){miu[i*pri[j]]=0;break;}
            miu[i*pri[j]]=-miu[i];
        }
        miu[i]+=miu[i-1];
    }
    scanf("%lld%lld",&L,&R);
    nn=L-1;memset(a,-1,sizeof(a));L=solve(L-1);
    nn=R;memset(a,-1,sizeof(a));R=solve(R);
    printf("%lld",R-L);
    return 0;
}

3.【bzoj3944】Sum

题意：求φ和μ的前缀和。

分析：无

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LL long long
using namespace std;
const int N=2e6;
LL T,n,nn,tot,pri[N+5];
LL miu[N+5],phi[N+5],mmiu[N+5],mphi[N+5];
bool np[N+5];
struct node{LL phi,miu;}ans,tmp;
node solve(LL n)
{
    if(n<=N)return (node){phi[n],miu[n]};
    if(~mphi[nn/n])return (node){mphi[nn/n],mmiu[nn/n]};
    LL ansp=1ll*n*(n+1)/2,ansm=1,pos;
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);tmp=solve(n/i);
        ansp=ansp-(pos-i+1)*tmp.phi;
        ansm=ansm-(pos-i+1)*tmp.miu;
    }
    mphi[nn/n]=ansp;mmiu[nn/n]=ansm;
    return (node){ansp,ansm};
}
void work()
{
    scanf("%lld",&n);nn=n;
    memset(mphi,-1,sizeof(mphi));
    ans=solve(n);
    printf("%lld %lld\n",ans.phi,ans.miu);
}
int main()
{
    phi[1]=miu[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i]){pri[++tot]=i;miu[i]=-1;phi[i]=i-1;}
        for(int j=1;i*pri[j]<=N;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0)
            {
                miu[i*pri[j]]=0;
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            miu[i*pri[j]]=-miu[i];
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
        miu[i]+=miu[i-1];
        phi[i]+=phi[i-1];
    }
    scanf("%lld",&T);
    while(T--)work();
    return 0;
}

4.【51nod1238】最小公倍数之和V3

题意：输出小于等于N的所有数，两两之间的最小公倍数之和。

分析：jiry_2の博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=5e6;
const int mod=1e9+7;
const int inv2=(mod+1)/2;
const int inv6=(mod+1)/6;
int tot,now,pri[N+5],mn[N+5];
bool np[N+5];
LL n,nn,sum,phi[N+5],a[N+5];
LL calc1(LL a)
{
    LL suma=a%mod;
    return suma*(suma+1)%mod*inv2%mod;
}
LL calc2(LL a,LL b)
{
    LL suma=(a-1)%mod;
    suma=suma*(suma+1)%mod*((2*suma+1)%mod)%mod*inv6%mod;
    LL sumb=b%mod;
    sumb=sumb*(sumb+1)%mod*((2*sumb+1)%mod)%mod*inv6%mod;
    return (sumb-suma+mod)%mod;
}
LL calc3(LL a,LL b)
{
    LL suma=(a-1)%mod;
    suma=suma*(suma+1)%mod*inv2%mod;
    suma=suma*suma%mod;
    LL sumb=b%mod;
    sumb=sumb*(sumb+1)%mod*inv2%mod;
    sumb=sumb*sumb%mod;
    return (sumb-suma+mod)%mod;
}
LL solve(LL n)
{
    if(n<=N)return phi[n];
    if(~a[nn/n])return a[nn/n];
    LL ans=0,pos;
    for(LL i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans+calc1(n/i)*calc3(i,pos))%mod;
    }
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans-solve(n/i)*calc2(i,pos)%mod+mod)%mod;
    }
    return a[nn/n]=ans;
}
int main()
{
    scanf("%lld",&n);
    nn=n;phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i])
        {
            pri[++tot]=i;
            phi[i]=(1ll*i*(i-1)%mod+1)%mod;
            mn[i]=i;
        }
        for(int j=1;i*pri[j]<=N;j++)
        {
            now=i*pri[j];
            np[now]=true;
            if(i%pri[j]==0)
            {
                mn[now]=mn[i]*pri[j];
                phi[now]=(phi[i]+phi[i/mn[i]]*mn[now]%mod*(mn[i]*(pri[j]-1)%mod))%mod;
                break;
            }
            phi[now]=phi[i]*(1ll*pri[j]*(pri[j]-1)%mod+1)%mod;
            mn[now]=pri[j];
        }
    }
    for(int i=1;i<=N;i++)phi[i]=phi[i]*i%mod;
    for(int i=2;i<=N;i++)phi[i]=(phi[i]+phi[i-1])%mod;
    memset(a,-1,sizeof(a));
    printf("%lld\n",solve(n));
    return 0;
}

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LL long long
using namespace std;
const int N=5e6;
const int mod=1e9+7;
const int inv2=(mod+1)/2;
const int inv6=(mod+1)/6;
int tot,pri[N+5];
LL n,nn;
LL phi[N+5],a[N+5];
bool np[N+5];
LL calc(LL a,LL b)
{
    LL suma=(a-1)%mod;
    suma=suma*(suma+1)%mod*((2*suma+1)%mod)%mod*inv6%mod;
    LL sumb=b%mod;
    sumb=sumb*(sumb+1)%mod*((2*sumb+1)%mod)%mod*inv6%mod;
    return (sumb-suma+mod)%mod;
}
LL solve(LL n)
{
    if(n<=N)return phi[n];
    if(~a[nn/n])return a[nn/n];
    LL ans=(n%mod)*((n+1)%mod)%mod*inv2%mod,pos;
    ans=ans*ans%mod;
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans-solve(n/i)*calc(i,pos)%mod+mod)%mod;
    }
    return a[nn/n]=ans;
}
int main()
{
    scanf("%lld",&n);
    nn=n;phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i]){pri[++tot]=i;phi[i]=i-1;}
        for(int j=1;i*pri[j]<=N;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
        phi[i]=phi[i]*i%mod*i%mod;
        phi[i]=(phi[i]+phi[i-1])%mod;
    }
    memset(a,-1,sizeof(a));
    LL pos,ans=0,sum;
    for(LL i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        sum=(solve(pos)-solve(i-1)+mod)%mod;
//        printf("   %lld %lld\n",pos,solve(pos));
        sum=sum*(((n/i)%mod)*((n/i+1)%mod)%mod*inv2%mod)%mod;
        ans=(ans+sum)%mod;
    }
    printf("%lld\n",ans);
//    for(int i=1;i<=n;i++)printf("%lld\n",phi[i]);
    return 0;
}

5.【51nod1237】最大公约数之和

题意：输出小于等于N的所有数，两两之间的最大公约数之和。

分析：无

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=5e6;
const int mod=1e9+7;
const int inv=(mod+1)/2;
int tot,pri[N+5];
LL n,nn,ans,pos,phi[N+5],a[N+5];
bool np[N+5];
LL calc(LL a,LL b)
{
    LL suma=(a-1)%mod;
    suma=suma*(suma+1)%mod*inv%mod;
    LL sumb=b%mod;
    sumb=sumb*(sumb+1)%mod*inv%mod;
    return (sumb-suma+mod)%mod;
}
LL solve(LL n)
{
    if(n<=N)return phi[n];
    if(~a[nn/n])return a[nn/n];
    LL ans=(n%mod)*((n+1)%mod)%mod*inv%mod,pos;
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans-(pos-i+1)%mod*solve(n/i)%mod+mod)%mod;
    }
    return a[nn/n]=ans;
}
int main()
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i]){pri[++tot]=i;phi[i]=i-1;}
        for(int j=1;i*pri[j]<=N;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
        phi[i]=(phi[i]+phi[i-1])%mod;
    }
    scanf("%lld",&n);nn=n;
    memset(a,-1,sizeof(a));
    for(LL i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans+calc(i,pos)*solve(n/i)%mod)%mod;
    }
    n%=mod;n=n*(n+1)%mod*inv%mod;
    ans=(ans*2%mod-n+mod)%mod;
    printf("%lld",ans);
    return 0;
}

6.【51nod1227】平均最小公倍数

题意：见原题

分析：无

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=1e6;
const int mod=1e9+7;
const int inv2=(mod+1)/2;
const int inv6=(mod+1)/6;
int tot,pri[N+5];
LL a,b,nn,pos,sum,phi[N+5],m[N+5];
bool np[N+5];
LL calc(LL a,LL b)
{
    LL suma=(a-1)%mod;
    suma=suma*(suma+1)%mod*inv2%mod;
    LL sumb=b%mod;
    sumb=sumb*(sumb+1)%mod*inv2%mod;
    return (sumb-suma+mod)%mod;
}
LL solve(LL n)
{
    if(n<=N)return phi[n];
    if(~m[nn/n])return m[nn/n];
    LL ans=n%mod,pos;
    ans=ans*(ans+1)%mod*((ans*2+1)%mod)%mod*inv6%mod;
    for(LL i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        sum=calc(i,pos)*solve(n/i)%mod;
        ans=(ans-sum+mod)%mod;
    }
    return m[nn/n]=ans;
}
LL work(LL n)
{
    LL ans=0;
    for(LL i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        sum=((pos-i+1)%mod)*solve(n/i)%mod;
        ans=(ans+sum)%mod;
//        printf("%lld\n",ans);
    }
    return (ans+n)%mod*inv2%mod;
}
int main()
{
    phi[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i]){pri[++tot]=i;phi[i]=i-1;}
        for(int j=1;i*pri[j]<=N;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];break;}
            phi[i*pri[j]]=phi[i]*(pri[j]-1);
        }
        phi[i]=phi[i]*i%mod;
        phi[i]=(phi[i]+phi[i-1])%mod;
    }
    scanf("%lld%lld",&a,&b);
    memset(m,-1,sizeof(m));nn=a;a=work(a-1);
    memset(m,-1,sizeof(m));nn=b;b=work(b);
    printf("%lld",(b-a+mod)%mod);
    return 0;
}

7.【bzoj4176】Lucas的数论

题意：见原题

分析：PoPoQQQの博客

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define LL long long
using namespace std;
const int N=6e6;
const int mod=1e9+7;
int S,tot,pri[N+5],miu[N+5],a[N+5];
int n,nn,pos,n2,pos2,sum,ans;
bool np[N+5];
int findmiu(int n)
{
    if(n<=S)return miu[n];
    if(~a[nn/n])return a[nn/n];
    int ans=1,pos;
    for(int i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans-(LL)(pos-i+1)*findmiu(n/i)%mod+mod)%mod;
    }
    return a[nn/n]=ans;
}
int solve(int a,int b)
{
    a=findmiu(a-1);b=findmiu(b);
    return (b-a+mod)%mod;
}
int main()
{
    scanf("%d",&n);
    S=ceil(pow(n,0.75)-1e-7)+1e-7;
    miu[1]=1;
    for(int i=2;i<=S;i++)
    {
        if(!np[i]){pri[++tot]=i;miu[i]=-1;}
        for(int j=1;i*pri[j]<=S;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0){miu[i*pri[j]]=0;break;}
            miu[i*pri[j]]=-miu[i];
        }
        miu[i]+=miu[i-1];
    }
    nn=n;memset(a,-1,sizeof(a));
    for(int i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);n2=n/i;sum=0;
        for(int j=1;j<=n2;j=pos2+1)
        {
            pos2=n2/(n2/j);
            sum=(sum+(LL)(pos2-j+1)*(n2/j)%mod)%mod;
        }
        sum=(LL)sum*sum%mod;
        ans=(ans+(LL)solve(i,pos)*sum%mod)%mod;
    }
    printf("%d",ans);
    return 0;
}

8.【51nod1220】约数之和

题意：见原题

分析：无。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long
using namespace std;
const int N=6e6;
const int mod=1e9+7;
int n,S,nn,tot,pos,ans;
int pri[N+5],miu[N+5],a[N+5];
bool np[N+5];
void pre()
{
    S=ceil(pow(n,0.75)-1e-7)+1e-7;
    miu[1]=1;
    for(int i=2;i<=S;i++)
    {
        if(!np[i]){pri[++tot]=i;miu[i]=-1;}
        for(int j=1;i*pri[j]<=S;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0){miu[i*pri[j]]=0;break;}
            miu[i*pri[j]]=-miu[i];
        }
        miu[i]=1ll*miu[i]*i%mod;
        miu[i]=(miu[i-1]+miu[i])%mod;
        miu[i]=(miu[i]+mod)%mod;
    }
}
int calc(int a,int b)
{
    a=1ll*a*(a+1)/2%mod;
    b=1ll*b*(b+1)/2%mod;
    return (b-a+mod)%mod;
}
int find(int n)
{
    if(n<=S)return miu[n];
    if(~a[nn/n])return a[nn/n];
    int ans=1,pos,sum;
    for(int i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        sum=1ll*calc(i-1,pos)*find(n/i)%mod;
        ans=(ans-sum+mod)%mod;
    }
    return a[nn/n]=ans;
}
int solve(int a,int b)
{
    a=find(a-1);b=find(b);
    return (b-a+mod)%mod;
}
int work(int n)
{
    int ans=0,pos;
    for(int i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans+1ll*calc(i-1,pos)*(n/i)%mod)%mod;
    }
    return 1ll*ans*ans%mod;
}
int main()
{
    scanf("%d",&n);
    pre();nn=n;memset(a,-1,sizeof(a));
    for(int i=1;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans+1ll*solve(i,pos)*work(n/i)%mod)%mod;
    }
    printf("%d",ans);
    return 0;
}

9.【bzoj3512】DZY Loves Math IV

题意：见原题

分析：permuiの博客

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
#define LL long long
using namespace std;
const int N=1e6;
const int mod=1e9+7;
const int inv=(mod+1)/2;
int n,m,ans,tot,now;
int pri[N+5],phi[N+5],mn[N+5],a[N+5];
bool np[N+5];
map<int,int>M[N+5];
int power(int a,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;b>>=1;
    }
    return ans;
}
int getsum(int n)
{
    if(n<=N)return phi[n];
    if(~a[m/n])return a[m/n];
    int ans=1ll*n*(n+1)%mod*inv%mod,pos;
    for(int i=2;i<=n;i=pos+1)
    {
        pos=n/(n/i);
        ans=(ans-1ll*(pos-i+1)*getsum(n/i)%mod+mod)%mod;
    }
    return a[m/n]=ans;
}
int Phi(int n)
{
    int ans=1,num;
    if(n<=N){ans=(phi[n]-phi[n-1]+mod)%mod;return ans;}
    int m=sqrt(n)+0.5;
    for(int i=1;pri[i]<=m;i++)
    {
        if(n%pri[i])continue;
        num=0;while(n%pri[i]==0)num++,n/=pri[i];
        ans*=power(pri[i],num-1)*(pri[i]-1);
    }
    if(n!=1)ans*=(n-1);
    return ans;
}
int S(int n,int m)
{
    if(m==0)return 0;
    if(n==1)return getsum(m);
    if(M[n][m])return M[n][m];
    int ans=0;
    for(int i=1;i*i<=n;i++)
    {
        if(n%i)continue;
        int j=n/i;ans=(ans+1ll*Phi(j)*S(i,m/i)%mod)%mod;
        if(i!=j)ans=(ans+1ll*Phi(i)*S(j,m/j)%mod)%mod;
    }
    return M[n][m]=ans;
}
int main()
{
    phi[1]=mn[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!np[i])pri[++tot]=i,phi[i]=i-1,mn[i]=i;
        for(int j=1;i*pri[j]<=N;j++)
        {
            now=i*pri[j];np[now]=true;
            if(i%pri[j]==0)
            {
                phi[now]=phi[i]*pri[j];
                mn[now]=mn[i];break;
            }
            phi[now]=phi[i]*(pri[j]-1);
            mn[now]=mn[i]*pri[j];
        }
        phi[i]=(phi[i-1]+phi[i])%mod;
    }
    memset(a,-1,sizeof(a));
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;i++)
        ans=(ans+1ll*i/mn[i]*S(mn[i],m)%mod)%mod;
    printf("%d",ans);
    return 0;
}

【洲阁筛】

〖相关资料〗

《洲阁筛学习》

〖相关题目〗

1.【spoj】DIVCNT3 - Counting Divisors (cube)

题意：见原题

分析：Monster_Yiの博客

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#define LL long long
using namespace std;
const int N=320000;
const int mx=316240;
const int p=1000003;
int tot,id,pri[N],num[N],ci[N],mn[N];
int cnt,first[p],w[p],nex[p],D[p];
LL T,n,ans1,f[N],f2[p],sum[N],to[p],d[p],g[p];
LL read()
{
    LL x=0,f=1;char c=getchar();
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
void ins(int x,LL y)
{
    int id=y%p;to[++cnt]=y;w[cnt]=x;
    nex[cnt]=first[id];first[id]=cnt;
}
int idx(LL x)
{
    for(int i=first[x%p];i;i=nex[i])
        if(to[i]==x)return w[i];
    return 0;
}
void work1()
{
    cnt=id=0;int k;
    for(LL i=1;i<=n;i=n/(n/i)+1)first[n/i%p]=0;
    for(LL i=1;i<=n;i=n/(n/i)+1)
        ins(++id,n/i),d[id]=g[id]=n/i,D[id]=0;
    for(int i=1;i<=tot;i++)
        for(int j=1;j<=id&&1ll*pri[i]*pri[i]<=d[j];j++)
            k=idx(d[j]/pri[i]),g[j]-=g[k]-(i-1-D[k]),D[j]=i;
}
void work2()
{
    for(int i=tot;i>=1;i--)
        for(int j=1;j<=id&&1ll*pri[i]*pri[i]<=d[j];j++)
        {
            if(pri[i+1]>d[j])f2[j]=1;
            else if(1ll*pri[i+1]*pri[i+1]>d[j])
                f2[j]=(num[min(1ll*mx,d[j])]-num[pri[i+1]-1])*4+1;
            for(LL pi=pri[i],c=1;pi<=d[j];pi*=pri[i],c++)
            {
                LL k=d[j]/pi,k2;
                if(pri[i+1]>k)k2=1;
                else if(1ll*pri[i+1]*pri[i+1]>k)
                    k2=(num[min(1ll*mx,k)]-num[pri[i+1]-1])*4+1;
                else k2=f2[idx(k)];
                f2[j]+=k2*(3*c+1);
            }
        }
}
LL solve(LL n)
{
    if(n<=mx)return sum[n];
    ans1=0;work1();work2();
    for(int i=1,pos;i<=mx;i=pos+1)
    {
         int j=idx(n/i);LL k;
         if(pri[tot+1]>n/i)k=0;
         else k=g[j]-(tot-D[j]+1);
         ans1+=(sum[pos=min(1ll*mx,n/(n/i))]-sum[i-1])*k*4;
    }
    return ans1+f2[1];
}
int main()
{
    f[1]=sum[1]=1;
    for(int i=2;i<=mx;i++)
    {
        num[i]=num[i-1]+(!f[i]);
        if(!f[i])pri[++tot]=i,f[i]=4,mn[i]=i,ci[i]=1;
        for(int j=1,now;(now=i*pri[j])<=mx;j++)
        {
            if(i%pri[j]==0)
            {
                ci[now]=ci[i]+1;
                f[now]=f[i/mn[i]]*(ci[now]*3+1);
                mn[now]=mn[i]*pri[j];
                break;
            }
            f[now]=f[i]*4;mn[now]=pri[j];ci[now]=1;
        }
        sum[i]=(sum[i-1]+f[i]);
    }
    pri[tot+1]=316241;
    T=read();
    while(T--)
    {
        n=read();
        printf("%lld\n",solve(n));
    }
    return 0;
}

2.【51nod 1184】第N个质数

题意：给出一个数N，求第N个质数。

分析：WerKeyTom_FTDの博客

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define LL long long 
using namespace std;
const int N=12500000;
int n,tot,f[N+5],pri[N+5];
LL l,r,mid,ans,tmp;
bool np[N+5];
LL g(LL n,int m)
{
    if(!m)return n;
    if(m==1)return n-n/2;
    if(n<=N)
    {
        if(m>=f[n])return 1;
        if(m>=f[(int)sqrt(n)])return f[n]-m+1;
    }
    return g(n,m-1)-g(n/pri[m],m-1);
}
LL calc(LL n)
{
    if(n<=N)return f[n];
    tmp=ceil(sqrt(n));
    return f[tmp]+g(n,f[tmp])-1;
}
int main()
{
    for(int i=2;i<=N;i++)
    {
        if(!np[i])pri[++tot]=i;
        for(int j=1;i*pri[j]<=N;j++)
        {
            np[i*pri[j]]=true;
            if(i%pri[j]==0)break;
        }
        f[i]=f[i-1]+(!np[i]);
    }
    scanf("%d",&n);
    l=1;r=22801763489;
    if(n<=500000000)r=11037271757;
    else l=11037271758;
    if(n<=750000000)r=16875026921;
    else l=16875026922;
    if(n<=850000000)r=19236701629;
    else l=19236701630;
    if(n<=900000000)r=20422213579;
    else l=20422213580;
    if(n<=940000000)r=22801763489;
    else l=22801763490;
    if(n<=970000000)r=22086742277;
    else l=22086742278;
    if(n<=990000000)r=22563323957;
    else l=22563323958;
    while(l<=r)
    {
        mid=(l+r)>>1;
        if(calc(mid)>=n)ans=mid,r=mid-1;
        else l=mid+1;
    }
    printf("%lld",ans);
    return 0;
}