浅谈 MIN_25 筛

在亚线性时间内求出一些积性函数前缀和的算法，这种递推还是太神了%%%。

一般适用条件

• 函数 \(f(x)\) 是积性函数。
• 函数 \(f(x)\) 是一个关于 \(x\) 可以快速求值的完全积性函数之和。（如多项式）
• 函数 \(f(p^k)\) 可以快速求值。（\(p\in prime\)）

算法过程

以 P5325 为例，讨论 MIN_25 筛的过程。

我们将质数和合数的贡献分类讨论。

设完全积性函数 \(f'\) 在质数处取值与 \(f\) 相同。

对于质数

设一个函数 \(g\)，满足 \(g(n,j) = \sum_{i=1}^{n} f'(i) \cdot [i \in \text{prime} \lor \min p|i > P_j]\)。（其实就是 \(1-n\) 的所有数被前 \(1\) 到 \(j\) 的质数筛去之后剩下数的和）

由于当 \(P_j\times P_j > n\) 之后没有筛去的贡献，于是考虑埃氏筛的过程，有转移方程：

\[g(n,j) = \begin{cases} g(n,j-1) & (P_j^2 > n) \\ g(n,j-1) - f'(P_j) \left[ g\left(\frac{n}{P_j}, j-1\right) - \sum_{i=1}^{j-1} f'(P_i) \right] & (P_j^2 \leq n) \end{cases} \]

相当于是将 \(P_j\) 的倍数都除去一个 \(P_j\)，考虑 \(k\times P_j\) 的 \(k\)。由于一些 \(k\times P_j\) 会在之前被 \(k\) 的因数筛掉，所以还要减掉重复计算的贡献，也就是后面的 \(\sum_{i=1}^{j-1} f'(P_i)\)。

这个时候，\(f'\) 作为完全积性函数的重要性就体现出来了，由于 \(f'(P_j) g\left(\frac{n}{P_j}, j-1\right)\) 本质上是将 \(f'(kP_j)\) 拆成了 \(f'(k)f'(P_j)\)，但是由于 \(k\) 与 \(P_j\) 不一定互质，所以 \(f'\) 是完全积性函数的条件是不可或缺的。对于此题，要将 \(f(x)\) 拆成两个单项式（都是完全积性函数）分别做，之后合并。

对于 \(g\) 的求法，发现第一维用到的只有 \(\frac{n}{x}\)，这个东西显然只有 \(2\sqrt{n}\) 个，对于 \(x\le\sqrt{n}\) 和 \(x>\sqrt{n}\) 分开拿数组存下来即可。

质数&合数

设函数 \(S\)，\(S(n,j) = \sum_{i=1}^{n} f(i) \cdot [\min p|i > P_j]\)。（与 \(g\) 的区别在于排除了前 \(1\) 到 \(j\) 的质数的贡献）

类似地，继续考虑埃式筛的过程，有：

\[S(n,j) = g(n,|P|) - \sum_{i=1}^{j} f'(P_i) + \sum_{k>j} \sum_{e=1}^{P_k^e \leq n} f(P_k^e) \left[ S\left( \frac{n}{P_k^e}, k \right) + (e>1) \right] \]

前一半是质数贡献，后一半是合数贡献。

后面 \((e>1)\) 的意义在于 \(S\left( \frac{n}{P_k^e}, k \right)\) 并不会把 \(1\) 算进去，所以要加上，但是当 \(e=1\) 的时候就不是合数了，所以当 \(e=1\) 的时候不能加一。

递归求解 \(S(n,0)\) 即为答案，记得加上 \(f(1)\)。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
inline ll read(){
	ll x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return x*f;
}
bool mbs;
const int maxn=2e5+20;
const int mod=1e9+7;
ll n,vis[maxn],prime[maxn],idx,sum1[maxn],sum2[maxn],g1[maxn],g2[maxn],id1[maxn],id2[maxn],tot,val[maxn];
inline void init(){
	for(int i=2;i<=100000;i++){
		if(!vis[i]) prime[++idx]=i;
		for(int j=1;j<=idx&&prime[j]*i<=100000;j++){
			vis[prime[j]*i]=1;
			if(i%prime[j]==0) break;
		}
	}
	for(int i=1;i<=idx;i++) sum1[i]=(sum1[i-1]+prime[i]*prime[i]%mod)%mod,sum2[i]=(sum2[i-1]+prime[i])%mod;
}
inline ll get(ll x){
	if(x<=(ll)sqrt(n)) return id1[x];
	return id2[n/x];
}
inline ll qpow(ll a,ll b){ll res=1;for(;b;b>>=1,a=a*a%mod) if(b&1) res=res*a%mod;return res;}
const ll inv2=qpow(2ll,mod-2);
const ll inv6=qpow(6ll,mod-2);
inline ll Sum2(ll x){x%=mod;return x*(x+1)%mod*inv2%mod;}
inline ll Sum1(ll x){x%=mod;return x*(x+1)%mod*(2*x%mod+1)%mod*inv6%mod;}
ll S(ll x,ll j){
	if(prime[j]>x) return 0;
	ll ans=(g1[get(x)]-g2[get(x)]+mod-(sum1[j]-sum2[j])+mod)%mod;
	for(int i=j+1;i<=idx&&prime[i]*prime[i]<=x;i++)  for(ll e=1,w=prime[i];w<=x;w*=prime[i],e++) 
		ans=(ans+w%mod*((w%mod-1)%mod)%mod*(S(x/w,i)+(e>1))%mod)%mod;
	return ans;
}
bool mbt;
int main(){
//	cerr<<(&mbs-&mbt)/1024.0/1024.0<<endl;
	init();
	n=read();
	for(ll l=1,r;l<=n;l=r+1){
		r=(n/(n/l)),val[++tot]=n/l;
		g1[tot]=(Sum1(val[tot])-1+mod)%mod,g2[tot]=(Sum2(val[tot])-1+mod)%mod;
		if(val[tot]<=(ll)sqrt(n)) id1[val[tot]]=tot;
		else id2[n/val[tot]]=tot;
	}
	for(int i=1;i<=idx;i++){
		for(int j=1;j<=tot&&prime[i]*prime[i]<=val[j];j++){
			g1[j]=(g1[j]-prime[i]*prime[i]%mod*(g1[get(val[j]/prime[i])]-sum1[i-1])%mod+mod)%mod;
			g2[j]=(g2[j]-prime[i]*(g2[get(val[j]/prime[i])]-sum2[i-1])%mod+mod)%mod;
		}
	}
	printf("%lld\n",(S(n,0)+1ll)%mod);
	return 0;
}

一些胡言乱语

可见 MIN_25 筛的精髓在于用类似递推的方式刻画埃式筛的过程。

至于为什么质数和合数要分开求，因为你筛合数的时候只需要 \(\le\sqrt{n}\) 的质数，如果将 \(e>1\) 直接改成 \(1\)，会导致 \(>\sqrt{n}\) 的质数没有计算。

一些应用

loj6235 区间素数个数

典型的应用，和上一题比较就是小巫见大巫了，容易有方程：

\[g(n,j)=g(n.j-1)-\left[g(\frac{n}{P_j},j-1)-(j-1)\right] \]

直接求就好。