min_25 筛学习小记

min_25筛

由 dalao min_25 发明的筛子，据说时间复杂度是极其优秀的 \(O(\frac {n^{\frac 3 4}} {\log n})\)，常数还小。

1. 质数 \(k\) 次方前缀和（基础）

求 \(\sum_{p \leq n}p^k\)

我们考虑一个 \(\rm DP\) 的思路：设 \(g(n,j)\) 为：

\[\sum_{i=1}^n[(\sum_{t=1}^j[p_t|i])=0] i^k \]

其实就是不大于 \(n\) 的，且不含有 \(p_1\) ~ \(p_j\) 中任意一个质因子的数的 \(k\) 次方之和。

我们考虑从 \(g(n,j-1)\) 转移上来。

发现 \(g(n,j)\) 相对于 \(g(n,j-1)\) 减少的就是含有 \(p_k\) 质因子的数的 \(k\) 次方之和，于是我们得到了转移方程：

\[g(n,k) = g(n,k-1)-p_j^kg(\lfloor \frac n {p_j} \rfloor,j) \]

但是，这个方程似乎有点儿问题？

我们把 \(p_j\) 以内的质数算重了，所以要加上。

最终的方程就是：

\[g(n,k) = g(n,k-1)-p_j^k(g(\lfloor \frac n {p_j} \rfloor,j)-g(p_{j-1},j-1) \]

不过如果我们要维护 \(O(n)\) 个 \(g\) 复杂度开销是很大的，所以考虑玄学优化去掉不重要的东西。

首先我们注意到，在 \(p_j > \sqrt n\) 是，对答案的贡献只有 \(p_j^k\)，所以我们可以先只考虑 \(p_j \leq \sqrt n\) 的情况。

\(g(p_{j-1},j-1)\) 显然可以在 \(O(\sqrt n)\) 的复杂度内线性筛预处理出来，就只需要考虑 \(g(\lfloor \frac n {p_j} \rfloor,j)\)

然后根据整除分块的结论，我们只需要处理 \(O(\sqrt n)\) 个 \(g\) 就行了。。。

然后复杂度就降下来了，我也不会证，不过看上去似乎是 \(O(\sqrt n\log\log n)\) 的（？）

但是 \(n\) 太大需要离散化。。。

所以我们依照杜教筛的离散化，把大于 \(\sqrt n\) 的映射到 \(1\) ~ \(\sqrt n\)。

给一下参考代码（这里只算了质数的 \(1\) 次方和 \(2\) 次方之和）：

const int St=1e6+5,mod=1e9+7;
inline int Add(const int&a,const int&b){
    return a+b>=mod?a+b-mod:a+b;
}
inline int Get(const ll&k){
	return k<=sqr?id1[k]:id2[n/k];
}
inline void seive(const int&n){
    int i,j,x;zhi[1]=1;
    for(i=2;i<=n;++i){
        if(!zhi[i]){
            pri[++top]=i;
            sum1[top]=Add(sum1[top-1],i);
            sum2[top]=Add(sum2[top-1],1ll*i*i%mod);
        }
        for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
            zhi[x]=1;
            if(!(i%pri[j]))break;
        }
    }
}
ll min_25(const ll&n){
    int i,j;
    seive(sqr=sqrt(n));
    for(ll L=1,R;L<=n;L=R+1){
        R=n/(n/L);
        w[++tot]=n/L;
        g1[tot]=w[tot]%mod;
        g2[tot]=g1[tot]*(g1[tot]+1)/2%mod*(g1[tot]*2+1)%mod*inv3%mod-1;
        g1[tot]=g1[tot]*(g1[tot]+1)/2%mod-1;
        if(n/L<=sqr)id1[n/L]=tot;
        else id2[n/(n/L)]=tot;
    }
    for(i=1;i<=top;++i){
        for(j=1;j<=tot&&pri[i]*pri[i]<=w[j];++j){
            ll k=Get(w[j]/pri[i]);
            g1[j]-=pri[i]*(g1[k]-sum1[i-1]+mod)%mod;
            g2[j]-=pri[i]*pri[i]%mod*(g2[k]-sum2[i-1]+mod)%mod;
            if(g1[j]<=0)g1[j]+=mod;
            if(g2[j]<=0)g2[j]+=mod;
        }
    }
}

在实现的时候将 \(g\) 的第二维滚掉了qwq

2.min_25筛

没错前面的都是铺垫（雾）

让我们看看我们要求的问题：

\[\sum_{i=1}^nf(i) \]

其中 \(f\) 是一个低阶多项式，其在质数乘方处的取值能够快速算出。

将答案分成两部分考虑：

\[\sum_{p \leq n}f(p)+\sum_{p^e \leq n,p \leq \sqrt n}f(p^e)(\sum_{i=1 \And minp > p}^{\lfloor \frac n {p^e} \rfloor}f(i)) \]

其中 \(minp\) 是 \(i\) 最小的质因子。

按照类似上面 \(\rm DP\) 的套路，设 \(S(n,x)\) 是不大于 \(n\) 的，且不含有 \(p_1\) ~ \(p_x\) 中任何一个因子的 \(f\) 之和。

显然答案是 \(S(n,0)\)。

再将答案分成两部分，一部分是大于 \(p_x\) 的质数的 \(f\) 之和，另一部分就是剩下的。

第一部分，也就是 \(g(n) - g(p_x)\)。

对照着上面的式子，我们可以得到以下的转移方程

\[S(n,x) = g(n) - g(p_x) +\sum_{p_k^e \leq n,k > x}f(p_k^e)(S(\lfloor \frac n {p_k^e} \rfloor,k)+[e!=1]) \]

然后还有一些细节，就是关于 \(1\) 的问题。。。

由于是照着 \(\rm wucstdio\) dalao 的题解学的，所以我也没有算 \(1\)，而是在结束时加上。

例题

板子

\[p^k(p^k-1) = p^{2k} - p^k \]

在质数处就是 \(p^2 - p\)

只需要筛出质数的 \(1\) 次方前缀和和 \(2\) 次方前缀和即可。

不要问我为什么上面的代码筛的也是 \(1\) 次和 \(2\) 次前缀和，我懒总行了吧

#include<cstdio>
#include<cmath>
typedef long long ll;
const int St=1e6+5,mod=1e9+7,inv3=333333336;
int sqr,top,tot,zhi[St],id1[St],id2[St];
ll n,w[St],g1[St],g2[St],pri[St],sum1[St],sum2[St];
inline int Add(const int&a,const int&b){
    return a+b>=mod?a+b-mod:a+b;
}
inline int Get(const ll&k){
	return k<=sqr?id1[k]:id2[n/k];
}
inline void seive(const int&n){
    int i,j,x;zhi[1]=1;
    for(i=2;i<=n;++i){
        if(!zhi[i]){
            pri[++top]=i;
            sum1[top]=Add(sum1[top-1],i);
            sum2[top]=Add(sum2[top-1],1ll*i*i%mod);
        }
        for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
            zhi[x]=1;
            if(!(i%pri[j]))break;
        }
    }
}
ll S(const ll&n,const int&k){
    if(pri[k]>=n)return 0;
    ll x=Get(n),ans=Add((g2[x]-g1[x]+sum1[k]-sum2[k])%mod,mod);
    for(int i=k+1;i<=top&&pri[i]*pri[i]<=n;++i){
        ll p=pri[i];
        for(int e=1;p<=n;++e,p*=pri[i]){
            ll id=p%mod;
            ans=Add(ans,id*(id-1)%mod*((e!=1)+S(n/p,i))%mod);
        }
    }
    return ans;
}
ll min_25(const ll&n){
    int i,j;
    seive(sqr=sqrt(n));
    for(ll L=1,R;L<=n;L=R+1){
        R=n/(n/L);
        w[++tot]=n/L;
        g1[tot]=w[tot]%mod;
        g2[tot]=g1[tot]*(g1[tot]+1)/2%mod*(g1[tot]*2+1)%mod*inv3%mod-1;
        g1[tot]=g1[tot]*(g1[tot]+1)/2%mod-1;
        if(n/L<=sqr)id1[n/L]=tot;
        else id2[n/(n/L)]=tot;
    }
    for(i=1;i<=top;++i){
        for(j=1;j<=tot&&pri[i]*pri[i]<=w[j];++j){
            ll k=Get(w[j]/pri[i]);
            g1[j]-=pri[i]*(g1[k]-sum1[i-1]+mod)%mod;
            g2[j]-=pri[i]*pri[i]%mod*(g2[k]-sum2[i-1]+mod)%mod;
            if(g1[j]<=0)g1[j]+=mod;
            if(g2[j]<=0)g2[j]+=mod;
        }
    }
    return Add(S(n,0),1);
}
signed main(){
    scanf("%lld",&n);
    printf("%d",min_25(n));
}

因为某些原因，对这个板子卡了常数，下面放出卡常后的代码：

#include<cstdio>
#include<cmath>
typedef unsigned ui;
typedef unsigned long long ull;
const ui M=1e5+5,mod=1e9+7;
ui S,top,pos[M],pri[M],p1[M],p2[M],s1[M],s2[M];ui g1[M<<1],g2[M<<1];double inv[M];ull n;
inline void Get(const ui&id,const ui&n){
	g1[id]=n*(n+1ull)/2%mod;g2[id]=333333336ull*g1[id]%mod*(2ull*n+1)%mod;--g1[id];--g2[id];
}
inline ui pow(ui a,ui b){
	ui ans(1);for(;b;b>>=1,a=1ull*a*a%mod)if(b&1)ans=1ull*ans*a%mod;return ans;
}
inline void sieve(){
	ui i,j,x;inv[1]=1;
	for(i=2;i<=S;++i){
		inv[i]=1./i+1e-15;
		if(!pos[i]){
			pri[pos[i]=++top]=i;
			s1[top]=(s1[top-1]+(p1[top]=i))%mod;
			s2[top]=(s2[top-1]+(p2[top]=1ull*i*i%mod))%mod;
		}
		for(j=1;j<=pos[i]&&(x=i*pri[j])<=S;++j)pos[x]=j;
	}
}
ui DFS(const ull&n,const ui&k,const ui&id){
	if(pri[k]>=n)return 0;
	ui e,x=n<=S?S+n:id,s;s=(mod+g2[x]-g1[x])%mod;s=(s+mod-(mod+s2[k]-s1[k])%mod)%mod;ull p;double I;
	for(ui i=k+1;i<=top&&pri[i]<=ui(n*inv[pri[i]]);++i){
		p=pri[i];s=(s+p*(p-1)%mod*DFS(ull(n*inv[p]),i,id*p))%mod;p*=pri[i];I=inv[pri[i]]*inv[pri[i]];
		for(e=2;p<=n;p*=pri[i],I*=inv[pri[i]],++e)x=p%mod,s=(s+x*(x-1ull)%mod*(1+DFS(n/p,i,id*p)))%mod;
	}
	return s;
}
ui min25(const ull&n){
	ui i,j,x,P1,P2,S1,S2;bool typ;S=sqrt(n);typ=1ull*S*S==n;sieve();
	for(i=1;i<S;++i)Get(i,ull(n*inv[i])%mod),Get(S+i,i);Get(S,ui(n*inv[S])%mod);if(!typ)Get(S<<1,S);
	for(i=1;i<=top;++i){
		const ui&P=pri[i];const double&X=n*inv[P],&INV=inv[P];P1=p1[i];P2=p2[i];S1=s1[i-1];S2=s2[i-1];
		for(j=1;(x=P*j)<=S;++j){
			g1[j]=(g1[j]+1ull*P1*(mod+S1-g1[x]))%mod;
			g2[j]=(g2[j]+1ull*P2*(mod+S2-g2[x]))%mod;
		}
		for(;1ull*P*(x=P*j)<=n&&j<=S;++j){
			x=X*inv[j];
			g1[j]=(g1[j]+1ull*P1*(mod+S1-g1[S+x]))%mod;
			g2[j]=(g2[j]+1ull*P2*(mod+S2-g2[S+x]))%mod;
		}
		for(j=S-typ;1ull*P*P<=j;--j){
			x=j*INV;
			g1[S+j]=(g1[S+j]+1ull*P1*(mod+S1-g1[S+x]))%mod;
			g2[S+j]=(g2[S+j]+1ull*P2*(mod+S2-g2[S+x]))%mod;
		}
	}
	return DFS(n,0,1)+1;
}
signed main(){
	scanf("%llu",&n);printf("%u",min25(n));
}

然后是一道很经典的三倍经验题DIVCNTK，做出这道题后 DIVCNT2 和 DIVCNT3 就是三倍经验。
设积性函数 \(f(p^e)=d((p^e)^k)=d(p^{ek})\)

在质数处的取值就是 $ k+1 $

#include<cstdio>
#include<cmath>
typedef unsigned long long ull;
const int M=2e6+5;
ull n,m,lim,sqr,tot,top,w[M],g[M],pri[M],zhi[M],id1[M],id2[M];
inline int Get(const ull&k){
	return k<=sqr?id1[k]:id2[n/k];
}
inline void sieve(const int&n){
	int i,j,x;zhi[1]=1;
	for(i=2;i<=n;++i){
		if(!zhi[i])pri[++top]=i;
		for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
			zhi[x]=1;
			if(!(i%pri[j]))break;
		}
	}
}
ull S(const ull&n,const int&k){
	if(pri[k]>=n)return 0;
	ull x=Get(n),ans=g[x]-k*(m+1);
	for(ull i=k+1;i<=lim&&pri[i]*pri[i]<=n;++i){
		for(ull e=1,p=pri[i];p<=n;++e,p*=pri[i]){
			ans+=(e*m+1)*(S(n/p,i)+(e!=1));
		}
	}
	return ans;
}
ull min_25(const ull&n){
	int i,j;
	sqr=sqrt(n);tot=0;lim=1;
	while(pri[lim]*lim[pri]<=n)++lim;--lim;
	for(ull L=1,R;L<=n;L=R+1){
		R=n/(n/L);
		w[++tot]=n/L;
		g[tot]=w[tot]-1;
		if(n/L<=sqr)id1[n/L]=tot;
		else id2[n/(n/L)]=tot;
	}
	for(i=1;i<=lim;++i){
		for(j=1;j<=tot&&pri[i]*pri[i]<=w[j];++j){
			ull k=Get(w[j]/pri[i]);
			g[j]-=g[k]-i+1;
		}
	}
	for(i=1;i<=tot;++i)g[i]*=(m+1);
	return S(n,0)+1;
}
signed main(){
	int i,T;
	sieve(2e6);
	scanf("%d",&T);
	for(i=1;i<=T;++i){
		scanf("%llu%llu",&n,&m);
		printf("%llu\n",min_25(n));
	}
}