【瞎口胡】数学 7 Min_25 筛

Min_25 筛用来求一类特殊数论函数的前缀和。

要求的函数 \(f\) 必须满足以下条件：

1. \(f\) 是积性函数。
2. 对于质数 \(p\)，\(f(p)\) 是一个关于 \(p\) 的简单低次多项式。
3. 对于质数 \(p\) 和正整数 \(e\)，\(f(p^e)\) 可以快速计算。

例题 Luogu P5325

转化题目

题目中的函数 \(f(p^k)=p^k(p^k-1)\) 显然满足要求。

即对于 \(x=p^k\)，\(f(x)=x^2-x\)，其中 \(p\) 是质数，\(k\) 是正整数。

显然我们可以将这样的函数拆成两个，即对于 \(x=p^k\)，把 \(f(x)\) 拆成 \(f_1(x)=x\) 和 \(f_2(x)=x^2\)，分别计算它们，求和后就是原来的 \(f\)。

在接下来的讨论中，我们设 \(f'(x)=x^k\)。分别将 \(k=1,2\) 带入并求和，就得到了答案。

观察到 \(f'\) 是完全积性函数，且在质数处和 \(f_k\) 取值相同。

来做一些奇怪的工作

设 \(\mathbb P\) 表示质数集， \(P_i\) 表示第 \(i\) 个质数。特别地，设 \(P_0=0\)。设 \(L(i)\) 表示 \(i\) 的最小质因子。考虑有这样一个数组 \(g\)，它满足以下定义：

\[g(n,j)=\sum \limits_{1 \leq i \leq n \and (i \in \mathbb P \or L(i)>P_j)} f'(i) \]

其中 \(j \leq Q\)，\(Q\) 是最大的满足 \(P_Q \leq \sqrt n\) 的正整数。

即，在 \(1\) 到 \(n\) 的所有 \(i\) 中，对所有满足 \(i\) 是质数或 \(i\) 的最小质因子大于 \(P_j\) 的 \(i\) 的 \(f'\) 值求和。

考虑 \(g(n,j)\) 的转移。如果 \(n<P_j^2\)，那么 \(n\) 的最小质因子必然不会为 \(P_j\)，所以这个限制不会对贡献产生影响，于是有 \(g(n,j)=g(n,j-1)\)。

如果 \(n \geq P_j^2\) 怎么办呢？这种情况下，\(g(n,j)\) 会在 \(g(n,j-1)\) 的基础上减少。观察到 \(P_{j-1},P_j\) 是连续的两个质数，于是从贡献中消失的 \(i\) 的最小质因子一定是 \(P_j\)。这一部分 \(i\) 的贡献即是 \(f'(P_j) \times( g(\left \lfloor \dfrac{n}{P_j} \right \rfloor,j-1)-g(P_{j-1},j-1))\)。为什么呢？观察式子开头的 \(f'(P_j)\)，这是因为 \(f'\) 是完全积性函数，且从贡献中消失的 \(i\) 的最小质因子一定是 \(P_j\)，我们可以将这一部分的贡献提出。\(g(\left \lfloor \dfrac{n}{P_j} \right \rfloor,j-1)\) 中的 \(i\)（姑且记作 \(i'\)）满足 \(L(i')>P_{j-1}\) 或 \(i' \in \mathbb P\)。

对于 \(L(i')>P_{j-1}\) 这一部分，\(i'\) 的最小质因子大于 \(P_{j-1}\)，那么 \(i' \times P_j\)（为什么要乘上 \(P_j\)？因为我们将它提到外面去了）的最小质因子只可能为 \(P_j\)，\(i= i' \times P_j\) 正是贡献中应该消失的 \(i\)。

为什么要加上 \(g(P_{j-1},j-1)\) 呢？这对应着第二部分 \(i' \in \mathbb P\)。对于质数 \(i'\)，仅当 \(i' \times P_j\) 的最小质因子为 \(P_j\) 时才被删掉。观察到对于小于 \(P_j\) 的质数 \(i'\)，我们不应该删去 \(i =i'\times P_j\) 的贡献。同时，观察到 \(i'\) 要小于等于 \(\left \lfloor \dfrac{n}{P_j} \right \rfloor\) 才会错误地删掉贡献。因为 \(P_j\) 一定不大于 \(\sqrt n\)，所以有 \(P_j \leq \left \lfloor \dfrac{n}{P_j} \right \rfloor\)，这个限制可以忽略。我们要加回来的就是小于 \(P_j\) 的质数的贡献，即前 \(j-1\) 个质数的 \(f'\) 值之和。由 \(g\) 的定义可知，这和 \(g(P_{j-1},j-1)\) 相等。

于是我们得到了 \(g\) 的递推式：

\[g(n,j) = \begin{cases} g(n,j-1) & P_j^2 > n \\ g(n,j-1)-f'(P_j) \times( g(\left \lfloor \dfrac{n}{P_j} \right \rfloor,j-1)-g(P_{j-1},j-1)) & P_j^2 \leq n\end{cases} \]

求解答案

设 \(S(n,x)\) 表示 \(\sum \limits_{1 \leq i \leq n} [L(i)> P_x] f_k(i)\)，即 \(1 \sim n\) 中所有最小质因子大于 \(P_x\) 的数的 \(f_k\) 值之和，注意是 \(f_k\) 而非 \(f'\)。

分类讨论：

• 如果符合条件的 \(i\) 是质数

  那么这一部分 \(i\) 之和为 \(g(n,Q)- \sum \limits_{i=1}^{x} f_k(P_i)\)。后面的和式可以预处理的时候做一个前缀和算出来。

• 如果符合条件的 \(i\) 是合数

  枚举最小质因子 \(p_j^e\)。这一部分的贡献即为

  \[\sum \limits_{p_j^e \leq n \and j > x }f_k(p^e) \times (S(\left \lfloor \dfrac{n}{p_j^e}\right \rfloor,j)+[e>1]) \]

  \(+[e>1]\) 是因为 \(p_j^e\) 在 \(e>1\) 的时候本身就是合数，应该被算进去。而 \(e=1\) 的时候就是质数，不在这个分类里面。

综合一下，我们可以得到

\[S(n,x) = g(n,Q)- \sum \limits_{i=1}^{x} f_k(P_i) + \sum \limits_{e \geq 1,j > x ,p_j^e \leq n}f_k(p_j^e) \times (S(\left \lfloor \dfrac{n}{p_j^e}\right \rfloor,j)+[e>1]) \]

这样问题就解决啦，答案就是 \(S(n,0)\)。值得一提的是，我们没有计算 \(f(1)\) 的值，因此还要加上 \(1\)。

时间复杂度据说是 \(O(\dfrac{n^{\frac34}}{\log n})\)，大概是 \(1000 \text{ms}\) 跑 \(n=10^{10}\)。

# include <bits/stdc++.h>

const int N=1000010,INF=0x3f3f3f3f,MOD=1e9+7,INV6=166666668,INV2=500000004;
typedef long long ll;

ll prime[N],psum;
bool vis[N];
ll g1[N],g2[N],prsum1[N],prsum2[N],w[N],wtot,idx1[N],idx2[N];

ll MAXN;

ll n;

inline int read(void){
	int res,f=1;
	char c;
	while((c=getchar())<'0'||c>'9')
		if(c=='-')f=-1;
	res=c-48;
	while((c=getchar())>='0'&&c<='9')
		res=res*10+c-48;
	return res*f;
}
inline void init(void){ // 预处理
	for(int i=2;i<=MAXN;++i){
		if(!vis[i]){
			prime[++psum]=i;
			prsum1[psum]=(prsum1[psum-1]+i)%MOD;
			prsum2[psum]=(prsum2[psum-1]+1ll*i*i)%MOD;
		}
		for(int j=1;i*prime[j]<=MAXN&&j<=psum;++j){
			vis[i*prime[j]]=true;
			if(i%prime[j]==0){
				break;
			}
		}
	}
	return;
}
inline ll getsum1(ll x){ // 记得取模
	x%=MOD;
	return x*(x+1)%MOD*INV2%MOD;
}
inline ll getsum2(ll x){
	x%=MOD;
	return x*(x+1)%MOD*(2ll*x+1)%MOD*INV6%MOD;
}
inline int getid(ll x){ // trick: n/x 只有 O(\sqrt n) 种取值，于是我们开大小为 2\sqrt n 的数组即可，idx1[x] 存储 x <= \sqrt n 时 n/x 的取值，idx2[n/x] 存储 x > \sqrt n 时 n/x 的取值
	return (x<=MAXN)?idx1[x]:idx2[n/x];
}
ll S(ll x,ll i){ // 递归计算 S
	if(prime[i]>=x){
		return 0;
	}
	int nowid=getid(x);
	ll res=(((g2[nowid]-g1[nowid])-(prsum2[i]-prsum1[i]))%MOD+MOD)%MOD;
	for(int j=i+1;j<=psum&&prime[j]*prime[j]<=x;++j){
		for(ll k=1,nowp=prime[j];nowp<=x;++k,nowp=nowp*prime[j]){
			ll truenowp=nowp%MOD;
			res=(res+(truenowp-1)*truenowp%MOD*(S(x/nowp,j)+(k>1)))%MOD;
		}
	}
	return res;
}
int main(void){
	scanf("%lld",&n);
	MAXN=1ll*sqrt(n);
	init();
	for(ll l=1,r;l<=n;l=r+1){
		ll val=n/l;
		r=n/val,w[++wtot]=val;
		g1[wtot]=getsum1(val)-1,g2[wtot]=getsum2(val)-1; // 预处理出 g(...,0)
		if(val<=MAXN){
			idx1[val]=wtot;
		}else{
			idx2[n/val]=wtot;
		}
	}
    // g1,g2 分别代表 k=1,2 时 g 数组的值
	for(int i=1;i<=psum;++i){ // 滚动计算 g
		for(int j=1;prime[i]*prime[i]<=w[j]&&j<=wtot;++j){
			int nowid=getid(w[j]/prime[i]);
			g1[j]=(g1[j]-prime[i]*(g1[nowid]-prsum1[i-1])%MOD+MOD)%MOD;
			g2[j]=(g2[j]-prime[i]*prime[i]%MOD*(g2[nowid]-prsum2[i-1])%MOD+MOD)%MOD;
		}
	}
	printf("%lld",(S(n,0)+1ll)%MOD);
	return 0;
}