自动搬运

查看原文

来自洛谷，原作者为

	渐变色 停止幻想

搬运于2025-08-24 22:45:16，当前版本为作者最后更新于2023-02-25 22:22:56，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

题目传送门

前言

尝试介绍下一种容易实现且常数很小的积性函数求和法。

Solution

式子的推导请见出题人 Leasier 的题解，进一步转化后得到原式 $=mS_{I_{m+1}}(n)$，其中 $I_n$ 表示 $I*I\,*···*\,I$ (卷 $n$ 次)，为积性函数。

对于积性函数 $f$，构造 $g$ 满足 $g(p^k)=f(p)^k$，则 $f=g*h$，易得 $h(p^k)=f(p^k)-f(p)f(p^{k-1})$。因此我们可以先算较容易求出的 $g$，再借助 $PN$ 得到 $f$。

回忆一下 min_25 筛的第一部分，$g$ 为满足 $g(p^k)=g(p)^k$ 的任一积性函数，令 $S(n, a)=\sum\limits_{i=1}^{n}[min_i>p_a \,or\;i\in P]\,g(i)$，有 :

$$S(n,a)=S(n,a-1)-g(p_a)(S(\frac{n}{p_a},a-1)-S(p_{a-1},a-1)) $$

那么其实它也能写成这样 :

$$S(n,a-1)=S(n,a)+g(p_a)(S(\frac{n}{p_a},a-1)-S(p_{a-1},a-1)) $$

也就是说对于我们要求的 $g$，我们只需先求出其在素数处的和（即得到 $S(m,\pi(\sqrt{n}))$，$m=\left\lfloor\dfrac{n}{i}\right\rfloor$，$i=1,2,···,\sqrt{n}$），便可通过后一式得到 $g$ 的块筛。

之后就可以在 $O(\sqrt n)$ 内得到 $S_f(n)$，或在 $O(\sqrt{n}\log^{1+\varepsilon} n)$ 内得到 $f$ 的块筛。

实现

远古代码直接粘贴过来了，欠佳。

#include <iostream>
#include <vector>
#include <cmath>
#include <functional>
#include <algorithm>
#include <time.h>
#include <fstream>
using namespace std;
using u32 = unsigned int;
using u64 = unsigned long long;
using i64 = long long;
double inv[100005];
u32 c[70][70];
vector<int> primes;
u32 solve(i64 N, u32 k) {
	int v = log2(N) + k;
	
	const auto div = [&](i64 n, int d) -> i64 {return n * inv[d]; };
	k++;
	v = sqrt(N + 0.5);
	for (int i = 1; i <= v; ++i) inv[i] = 1.0 / i * (1 + 1e-15);
	vector<u32> s(v + 1), l(v + 1);
	for (int i = 1; i <= v; ++i) s[i] = i - 1, l[i] = div(N, i) - 1;
	for (int p = 2; p <= v; ++p)
		if (s[p] != s[p - 1]) {
			primes.push_back(p);
			i64 q = i64(p) * p, M = div(N, p);
			u32 t0 = s[p - 1];
			int t = v / p, u = min<i64>(v, div(M, p));
			for (int i = 1; i <= t; ++i) l[i] -= l[i * p] - t0;
			for (int i = t + 1; i <= u; ++i) l[i] -= s[div(M, i)] - t0;
			for (int i = v, j = t; j >= p; --j)
				for (int l = j * p; i >= l; --i)
					s[i] -= s[j] - t0;
		}
	for (int i = 1; i <= v; ++i) s[i] *= k, l[i] *= k;
	for (auto it = primes.rbegin(); it != primes.rend(); ++it) {
		int p = *it;
		i64 q = i64(p) * p, M = div(N, p);
		int t = v / p, u = min<i64>(v, div(M, p));
		u32 t0 = s[p - 1];
		if(q <= v) {
			u32 s0 = k * k;
			for (int j = p, i = q; j < t; s0 = (s[++j] - t0) * k)
				for (int l = j * p + p; i < l; ++i) s[i] += s0;
			for (int i = t * p; i <= v; ++i) s[i] += s0;
		}
		for (int i = u; i > t; --i) l[i] += (s[div(M, i)] - t0) * k;
		for (int i = t; i >= 1; --i) l[i] += (l[i * p] - t0) * k;
	}
	for (int i = 1; i <= v; ++i) s[i] += 1, l[i] += 1;
	
	k--;
	const auto divide = [](i64 n, i64 d) -> i64 {return double(n) / d; };
	function< u32(i64, int, u32)> rec = [&] (i64 n, size_t beg, u32 coeff) -> u32 {
		if (!coeff) return 0;
		u32 ret = coeff * (n > v ? l[divide(N, n)] : s[n]);
		for (size_t i = beg; i < primes.size(); ++i) {
			int p = primes[i];
			i64 q = i64(p) * p;
			if (q > n) break;
			i64 nn = divide(n, q);
			for (u32 e = 1; nn; nn = div(nn, p), ++e)
				ret += rec(nn, i + 1, -coeff * e * c[e + k][e + 1]);
		}
		return ret;
	};
	return k * rec(N, 0, 1);
}
int main() {
	i64 N;
	u32 k;
	cin >> N >> k;
	cout << solve(N, k);
	return 0;
}