min25筛

min_25 筛

非常棒的博客

已知一积性函数 (f(x))，(f(p)) 为关于 (p) 的多项式，且 (f(p^k)) 可以快速算出来，求 (sum_{i=1}^nf(i))。

例：求 (sum_{i=1}^n sigma_0(i^k))，其中 (n = 10^{10})

step1 : 求质数的积性函数前缀和

求 (1 ... n) 以内的质数个数。

我们设 (f(n,m)) 表示 (1...n) 中的所有质数或者不含最小质因子为 (p_1...p_m) （(p_i) 表示第 (i) 个质数）的数的个数。

那么有：

[f(n,m) = f(n,m-1) + (f(lfloor frac{n}{p_m} floor,m-1) - (m-1)) ]

意义为：考虑从 (f(n,m-1)) 处转移过来，那么我们只需要筛掉最小质因子恰好是 (p_m) 的数。那么我们把所有要筛的数提出来个 (p_m)，剩下的数的最小质因子为 (p_1...p_{m-1})。注意删掉“质数"部分。

更普遍的情况。

类似地，我们设 (f(n,m)) 表示 (1...n) 中的所有质数或者不含最小质因子为 (p_1...p_m) （(p_i) 表示第 (i) 个质数）的数的函数和。这时候我们需要认为合数也有一个假的和质数一样的函数和。那么有：

[f(n,m) = f(n,m-1) + F(p_m)(f(lfloor frac{n}{p_m} floor,m-1) - psum(m-1)) ]

（函数重名了，用 (F(n)) 代替了积性函数；(psum) 为预处理出来的质数函数前缀和）

由于是多项式，我们可以拆成 (F(p) = p^k) 的形式，最终拼起来。这样，假函数成为完全积性函数，直接提出来 (F(p_m)) 就是对的了。

复杂度证明

见那篇博客。

发现 (p_m^2 le n)

可以积分证明，复杂度为 (O(n^{0.75} / log n))

step2 : 求所有数的积性函数和

我们设 (g(n,m)) 表示 (1...n) 中最小质因子大于 (p_m) 的数的（真）积性函数和。

那么有：

[g(n,m) = f(n,m) - psum(m-1) + sum_{p_k^e} F(p_k^e)(g(lfloor frac{n}{p^e} floor,k) + [e ot= 1]) ]

意义为：分质数，合数考虑贡献。质数的贡献我们已经搞好了，为 (f(n,m) - psum(m-1))，合数的话我们需要枚举最小质因子及其次数，即枚举 (k > m) 的 (p_k)，再枚举次数 (e)。由于是积性函数，我们可以直接将 (F(p_k^e)) 提出来。值得注意的是，形如 (p^k,k > 1) 的数我们既没有在质数处枚举到，也没有在合数处提出来 (p_k^e) 后枚举到，所以需要加上个 ([e ot= 1])。

复杂度：不会证。状态数和上面一样是 (O(n^{0.75} / log n)) 的，但是转移复杂度比较玄学，据说总复杂度是 (O(n^{0.75} / log n)) 的，也有的说是 (O(n^{1 - epsilon})) 的。

代码：

//divcntk
int block, itot, pos[N];
ull id[N], n, K, f[N];
inline int get_id(ull x) {
	return pos[x <= block ? x : block + n / x];
}
inline void get_f() {
	for (ull l = 1, r; l <= n; l = r + 1) {
		ull zhi = n / l;
		r = n / zhi;
		id[++itot] = zhi;
		pos[zhi <= block ? zhi : block + n / zhi] = itot;
		f[itot] = zhi - 1;
	}
	for (int j = 1; 1ull * pri[j] * pri[j] <= n; ++j) {
		for (int i = 1; id[i] >= 1ull * pri[j] * pri[j]; ++i) {
			f[i] = f[i] - f[get_id(id[i] / pri[j])] + j - 1;
		}
	}
}

ull get_g(ull n, int j) {
	if (pri[j] >= n)	return 0;
	ull res = (f[get_id(n)] - j) * (K + 1);//bug
	for (int k = j + 1; 1ull * pri[k] * pri[k] <= n; ++k) {
		ull tmp = pri[k];
		for (int e = 1; tmp <= n; ++e, tmp = tmp * pri[k]) {
			res += (K * e + 1) * (get_g(n / tmp, k) + (e > 1));
		}
	}
	return res;
}

inline void work() {
	read(n), read(K); block = sqrt(n);//bug
	get_f();
	ull res = get_g(n, 0) + 1;
	printf("%llu
", res);
}

min_25 筛