min25筛学习笔记

min25筛简介：用来求积性函数F(x)前缀和的，复杂度O(n/logn)，大概能求n<=10。

记一个数x的最小质因子为R(x)，所以当x不为质数时，R(x)<=√x这是废话。

首先求所有质数的F(x)和，下设g(i,j)=ΣF(x)，其中2<=x<=i，且x为质数或R(x)>pri[j]，其中pri[j]为第j个质数。其实，j的取值至多√n个显而易见，下面可以发现最终状态是g(i,tot)，其中tot为√n以内的质数个数。初始化g(i,0)，即将所有数视为质数处理，然后将合数筛去，剩下的即为所要求的值。

然后g的递推式也是显而易见，用g(i,j-1)推g(i,j)，可以得到g(i,j)=g(i,j-1)-F(pri[j])*(g(i/pri[j],j-1)-sF(j-1))，其中sF(x)=ΣF(pri[i])，其中1<=i<=x。

注：此时i/pri[j]可能会小于pri[j]，但此时显然g(i,j)=g(i,j-1)

第一维的大小只有2√n，显然应该离散存储。

然后下面考虑求所有数的F(x)前缀和，此时类似g的定义设f(i,j)=ΣF(x)，其中2<=x<=i，且x为质数或R(x)>=pri[j]。对于f(i,j)显然质数部分为g(i,tot)-sF(j-1)，合数部分，枚举最小因子及其次幂p使p>=pri[j]，然后F函数实际上是F(p)*f(i/p,p的序号+1)，不知道为啥，不需要记忆化。

例题：LOJ6053

#include<bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N=1e5+,mod=1e9+,inv2=;

ll n,P,cnt,id1[N],id2[N],w[*N],m,g[*N],h[*N],f[*N],pri[N],pre[N],vis[N];

void init()

{

    for(int i=;i<=1e5;++i)

    {

        if(!vis[i])pri[++cnt]=i,pre[cnt]=(pre[cnt-]+i)%mod;

        for(int j=;j<=cnt&&i*pri[j]<=1e5;++j)

        {

            vis[i*pri[j]]=;

            if(i%pri[j]==)break;

        }

    }

}

int getid(ll x){return x<=P?id1[x]:id2[n/x];}

void calc()

{

    ll v,r;

    for(ll l=;l<=n;l=r+)

    {

        v=n/l,r=n/v;

        if(v<=P)id1[v]=++m;else id2[r]=++m;

        w[m]=v;

        ll z=v%mod;

        g[m]=(z+)*(z-)%mod*inv2%mod;

        h[m]=z-;

    }

    for(ll i=;i<=cnt;++i)

    for(ll j=;j<=m&&pri[i]*pri[i]<=w[j];++j)

    {

        int op=getid(w[j]/pri[i]);

        g[j]=(g[j]-pri[i]*(g[op]-pre[i-])%mod)%mod;

        h[j]=(h[j]-(h[op]-i+))%mod;

    }

    for(int i=;i<=m;++i)g[i]=(g[i]+mod)%mod,f[i]=(f[i]+mod)%mod,f[i]=(g[i]-h[i])%mod;

}

ll ask(ll x,ll i)

{

    if(x<=||x<pri[i])return ;

    ll ans=f[getid(x)]-(pre[i-]-i+);

    if(i==)ans+=;

    for(int j=i;j<=cnt&&pri[j]*pri[j]<=x;++j)

    {

        ll r=pri[j];

        for(ll e=;r*pri[j]<=x;++e,r*=pri[j])

        ans=(ans+(pri[j]^e)*ask(x/r,j+)+(pri[j]^(e+)))%mod;

    }

    return ans;

}

int main()

{

    cin>>n;P=sqrt(n);

    init();

    calc();

    printf("%lld",((ask(n,)+)%mod+mod)%mod);

}

SPOJ34096 DIVCNTK

这题也太鬼畜了，1可以数论分块，2可以moblus反演，3可以洲阁筛。

#include<bits/stdc++.h>

using namespace std;

typedef unsigned long long ll;

const int N=1e5+;

bitset<N>vis;

ll n,k,lim,w[N<<],sp[N],g[N<<],h[N<<];

int m,num,tot,pri[N],id1[N],id2[N];

ll S(ll x,int y)

{

    if(x<=||pri[y]>x)return ;

    int t=x<=lim?id1[x]:id2[n/x];

    ll ret=g[t]+h[t]-(k+)*(y-);

    for(int i=y;i<=tot&&1ll*pri[i]*pri[i]<=x;i++)

    {

        ll p=pri[i];

        for(int j=;p*pri[i]<=x;p*=pri[i],j++)

        {

            t=x/p<=lim?id1[x/p]:id2[n/(x/p)];

            ret+=S(x/p,i+)*(k*j+)+k*(j+)+;

        }

    }

    return ret;

}

int main()

{

    for(int i=;i<N;i++)

    {

        if(!vis[i])pri[++tot]=i;

        for(int j=;j<=tot&&i*pri[j]<N;j++)

        {

            vis[i*pri[j]]=;

            if(i%pri[j]==)break;

        }

    }

    num=tot;

    int T;scanf("%d",&T);

    while(T--)

    {

        scanf("%lld%lld",&n,&k),lim=sqrt(n);

        m=;

        tot=upper_bound(pri+,pri+num+,lim)-pri-;

        for(ll i=,j;i<=n;i=j+)

        {

            j=n/(n/i),w[++m]=n/i;

            if(w[m]<=lim)id1[w[m]]=m;else id2[n/w[m]]=m;

            g[m]=(w[m]-)*k,h[m]=w[m]-;

        }

        for(int j=;j<=tot;j++)

        for(int i=;1ll*pri[j]*pri[j]<=w[i];i++)

        {

            int t=w[i]/pri[j]<=lim?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];

            g[i]-=g[t]-(j-)*k,h[i]-=h[t]-(j-);

        }

        printf("%llu\n",S(n,)+);

    }

}