min25筛学习笔记

min25筛简介：用来求积性函数F(x)前缀和的，复杂度O(n^0.75/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^k使p>=pri[j]，然后F函数实际上是F(p^k)*f(i/p^k,p的序号+1)，不知道为啥，不需要记忆化。

例题：LOJ6053

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e5+1000,mod=1e9+7,inv2=500000004;
ll n,P,cnt,id1[N],id2[N],w[2*N],m,g[2*N],h[2*N],f[2*N],pri[N],pre[N],vis[N];
void init()
{
    for(int i=2;i<=1e5;++i)
    {
        if(!vis[i])pri[++cnt]=i,pre[cnt]=(pre[cnt-1]+i)%mod;
        for(int j=1;j<=cnt&&i*pri[j]<=1e5;++j)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0)break;
        }
    }
}
int getid(ll x){return x<=P?id1[x]:id2[n/x];}
void calc()
{
    ll v,r;
    for(ll l=1;l<=n;l=r+1)
    {
        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+2)*(z-1)%mod*inv2%mod;
        h[m]=z-1;
    }
    for(ll i=1;i<=cnt;++i)
    for(ll j=1;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-1])%mod)%mod;
        h[j]=(h[j]-(h[op]-i+1))%mod;
    }
    for(int i=1;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<=1||x<pri[i])return 0;
    ll ans=f[getid(x)]-(pre[i-1]-i+1);
    if(i==1)ans+=2;
    for(int j=i;j<=cnt&&pri[j]*pri[j]<=x;++j)
    {
        ll r=pri[j];
        for(ll e=1;r*pri[j]<=x;++e,r*=pri[j])
        ans=(ans+(pri[j]^e)*ask(x/r,j+1)+(pri[j]^(e+1)))%mod;
    }
    return ans;
}
int main()
{
    cin>>n;P=sqrt(n);
    init();
    calc();
    printf("%lld",((ask(n,1)+1)%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+7;
bitset<N>vis;
ll n,k,lim,w[N<<1],sp[N],g[N<<1],h[N<<1];
int m,num,tot,pri[N],id1[N],id2[N];
ll S(ll x,int y)
{
    if(x<=1||pri[y]>x)return 0;
    int t=x<=lim?id1[x]:id2[n/x];
    ll ret=g[t]+h[t]-(k+1)*(y-1);
    for(int i=y;i<=tot&&1ll*pri[i]*pri[i]<=x;i++)
    {
        ll p=pri[i];
        for(int j=1;p*pri[i]<=x;p*=pri[i],j++)
        {
            t=x/p<=lim?id1[x/p]:id2[n/(x/p)];
            ret+=S(x/p,i+1)*(k*j+1)+k*(j+1)+1;
        }
    }
    return ret;
}
int main()
{
    for(int i=2;i<N;i++)
    {
        if(!vis[i])pri[++tot]=i;
        for(int j=1;j<=tot&&i*pri[j]<N;j++)
        {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0)break;
        }
    }
    num=tot;
    int T;scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%lld",&n,&k),lim=sqrt(n);
        m=0;
        tot=upper_bound(pri+1,pri+num+1,lim)-pri-1;
        for(ll i=1,j;i<=n;i=j+1)
        {
            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]-1)*k,h[m]=w[m]-1;
        }
        for(int j=1;j<=tot;j++)
        for(int i=1;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-1)*k,h[i]-=h[t]-(j-1);
        }
        printf("%llu\n",S(n,1)+1);
    }
}