【瞎口胡】数学 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) imes( g(left lfloor dfrac{n}{P_j} ight floor,j-1)-g(P_{j-1},j-1)))。为什么呢？观察式子开头的 (f'(P_j))，这是因为 (f') 是完全积性函数，且从贡献中消失的 (i) 的最小质因子一定是 (P_j)，我们可以将这一部分的贡献提出。(g(left lfloor dfrac{n}{P_j} ight floor,j-1)) 中的 (i)（姑且记作 (i')）满足 (L(i')>P_{j-1}) 或 (i' in mathbb P)。

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

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

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

[g(n,j) = egin{cases} g(n,j-1) & P_j^2 > n \ g(n,j-1)-f'(P_j) imes( g(left lfloor dfrac{n}{P_j} ight floor,j-1)-g(P_{j-1},j-1)) & P_j^2 leq nend{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 > i }f_k(p^e) imes (S(left lfloor dfrac{n}{p_j^e} ight floor,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_{p_j^e leq n and j > i }f_k(p^e) imes (S(left lfloor dfrac{n}{p_j^e} ight floor,j)+[e>1]) ]

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

时间复杂度据说是 (O(dfrac{n^{frac34}}{log n}))，大概是 (1000 ext{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) 种取值，于是我们开大小为 2sqrt 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;
}