BZOJ2627 JZPKIL

nb 莫比乌斯反演题，和常规套路不同。

题意：

$$\sum_{i = 1} ^ n \gcd(i, n) ^ x \text{lcm}(i, n) ^ y \bmod 10 ^ 9 + 7$$

$T(T\leq 100)$ 次询问，每次给出 $n, x, y$，$n\leq 10 ^ {18}$，$x, y\leq 3000$。

看到 $n\leq 10 ^ {18}$ 似乎不对劲，但还是先套路莫反。

$$\begin{aligned} &\sum_{i = 1} ^ n \sum_{d = 1} ^ n [\gcd(i, n) = d] d ^ x (\dfrac{in}d) ^ y\\ =& \sum_{d | n} \sum_{i = 1} ^ {\frac nd} d ^ x (in) ^ y [\gcd(i, \dfrac nd) = 1]\\ =& \sum_{d | n} d ^ x \sum_{i = 1} ^ n d ^ x (in) ^ y \sum_{s | i, s | \frac nd} \mu(s)\\ =& n ^ y \sum_{d | n} d ^ x \sum_{s | \frac nd}\mu(s) s ^ y \sum_{i = 1} ^ {\frac n{ds}} i ^ y \end{aligned}$$

做到这里，似乎再化开，令 $T = ds$ 没有意义了，因为至少需要枚举 $T$，少说也得整除分块，复杂度 $O(\sqrt n)$ 不可接受。

发现一个绊脚石是 $\sum_{i = 1} ^ {\frac n{ds}}i ^ y$，这个自然数幂导致我们没法快速计算，因为前面的函数 $d ^ x, \mu(s) s ^ y$ 都可以狄利克雷卷积，但这个没法。

考虑使用伯努利数，把这样的一个式子写成 $\dfrac n{ds}$ 的次幂和：

$$\begin{aligned} &\sum_{i = 1} ^ {\frac n {ds}} i ^ y\\ =& \sum_{k = 0} ^ y \binom{y + 1}k B_k ^ + (\dfrac n{ds}) ^ {y + 1 - k} \end{aligned}$$

其中 $B_k^+$ 和 $B_k$ 的伯努利数略有区别，$B_k^+$ 在 $B_1^+$ 时取的 $\dfrac 12$，而 $B_k$ 在 $B_1$ 取的是 $-\dfrac 12$。伯努利数可以在 $O(y ^ 2)$ 的时间内递推出前 $y$ 项，不再赘述，可以接受。

带入原式，大力计算：

$$\begin{aligned} & n ^ y \sum_{d | n} d ^ x \sum_{s | \frac nd}\mu(s) s ^ y \sum_{i = 1} ^ {\frac n{ds}} i ^ y\\ =& \dfrac{n ^ y}{y + 1} \sum_{d | n} d ^ x \sum_{s | \frac nd}\mu(s) s ^ y \sum_{k = 0} ^ y \binom{y + 1}k B_k^+ (\dfrac n{ds}) ^ {y + 1 - k}\\ =& \dfrac{n ^ y}{y + 1} \sum_{k = 0} ^ y \binom{y + 1}k B_k ^ + \sum_{d | n} d ^ x \sum_{s | \frac nd} \mu(s) s ^y (\dfrac{n}{ds}) ^ {y + 1 - k} \end{aligned}$$

枚举 $k$ 是可以接受的，那么现在考虑怎么计算后面的 $\sum_{d | n} d ^ x \sum_{e | \frac nd} \mu(s)s ^ y(\dfrac n{ds}) ^ {y + 1 - k}$。发现这个其实一个积性函数，是 $id ^ x \times (\mu\cdot id ^ y)\times id ^ {y + 1 - k}$，可以拆开每一个质数计算答案。而每一个质数计算的时候，相当于计算 $f(p ^ c)$，由于 $\mu$ 的有效位置只有 $1, p$ 两个位置，只需要考虑剩下的两个分配次幂即可。得到质因数的办法显然 Pollard-Rho，时间复杂度 $O(\sqrt p)$。

那么预处理时间复杂度为 $O(y ^ 2)$，单次询问时间复杂度 $O(n ^ {\frac 14} + y\cdot \text{poly}\log n)$，可以通过。

代码比较冗长，主要是 Pollard-Rho 太长了，就放一个求解函数。

void work()
{
	LL n;
	int x, y;
	std::cin >> n >> x >> y;
	std::vector<LL> tmpfac = Pollard_Rho(n);
	LL tmp = n;
	std::vector<PLI> fac;
	for (LL &x : tmpfac)
	{
		int cnt = 0;
		while (tmp % x == 0) cnt ++, tmp /= x;
		fac.push_back({x, cnt});
	}
	int res = 0;
	for (int k = 0; k <= y; ++ k)
	{
		int cur = 1;
		for (auto [p, t] : fac)
		{
			p %= Mod;
			int mul = 0;
			for (int i = 0; i <= t; ++ i)
				adj(mul += qpow(p, i * (y + 1 - k) + (t - i) * x) - Mod);
			for (int i = 0; i < t; ++ i)
				adj(mul -= qpow(p, i * (y + 1 - k) + y + (t - i - 1) * x));
			cur = (LL) cur * mul % Mod;
		}
		res = (res + (LL) cur * C(y + 1, k) % Mod * B[k]) % Mod;
	}
	res = (LL) res * qpow(n % Mod, y) % Mod * qpow(y + 1) % Mod;
	std::cout << res << '\n';
}