LOJ# 572. 「LibreOJ Round #11」Misaka Network 与求和（min25筛，杜教筛，莫比乌斯反演）

题意

求

$$\sum_{i = 1}^{n} \sum_{i = 1}^{n} f(\gcd(i, j))^k \pmod {2^{32}}$$

其中 $f(x)$ 为 $x$ 的次大质因子，重复的质因子计算多次。

特别的，定义 $f(1) = 0, f(p) = 0$ ，此处 $p$ 为质数。

题解

首先先莫比乌斯反演前几步。

$$ans = \sum_{d = 1}^{n} f(d)^k \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(x) (\lfloor \frac{n}{dx} \rfloor)^2$$

令 $T = dx$ 那么就化为

$$= \sum_{T = 1}^{n} (\lfloor \frac{n}{T} \rfloor)^2 \sum_{d | T} f(d)^k \mu(\frac{T}{d})$$

令 $f(d)^k = F(d)$ 那么就变成

$$= \sum_{T = 1}^{n} (\lfloor \frac{n}{T} \rfloor)^2 (F * \mu)(T)$$

那么我们整除分块后，只需要快速求 $(F * \mu)$ 的前缀和即可，令 $\displaystyle S(n) = \sum_{i = 1}^n (F * \mu)(i)$ 。

我们知道 $\mu * 1 = \epsilon$ ，由于狄里克雷卷积满足结合律，就有 $F * \mu * 1 = F$ 。

所以我们套上杜教筛的式子，就可以得到

$$S(n) = \sum_{i = 1}^n F(i) - \sum_{i = 2}^n S(\lfloor \frac n i \rfloor)$$

杜教筛好像不好算 $\displaystyle \sum_{i = 1}^n F(i)$ ，其实这是道 原题 。

min25 筛合数的部分 十分契合次大质因子的过程。

每次我们枚举了最小质因子的值，下一次递归计算 $S(n, i)$ 的时候，那么在 $P_i \sim P_{|P|}$ 的所有质数上次枚举过来的次大质因子一定都是 $P_{i - 1}$ 。

那么直接在此处加上贡献即可，不要忘记 $p^e$ 的贡献要加上去。

复杂度？？？如果考虑预处理前 $n^{2 / 3}$ 的答案，并且对于合数部分记忆化的话，应该比较正确。。。

但似乎直接实现虽然有点慢，但也能过？？？

代码

#include <bits/stdc++.h>

#define For(i, l, r) for (register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for (register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Rep(i, r) for (register int i = (0), i##end = (int)(r); i < i##end; ++i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl

using namespace std;

using ll = long long;
using ui = unsigned int;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

inline ui read() {
	ui x(0), sgn(1); char ch(getchar());
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') sgn = -1;
	for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
	return x * sgn;
}

void File() {
#ifdef zjp_shadow
	freopen ("572.in", "r", stdin);
	freopen ("572.out", "w", stdout);
#endif
}

const int N = 1e5 + 1e3;

int prime[N], pcnt; bitset<N> is_prime; 

ui Pow[N], k;

ui fpm(ui x, ui power) {
	ui res = 1;
	for (; power; power >>= 1, x *= x)
		if (power & 1) res *= x;
	return res;
}

void Linear_Sieve(int maxn) {
	is_prime.set();
	For (i, 2, maxn) {
		if (is_prime[i]) 
			prime[++ pcnt] = i, Pow[pcnt] = fpm(i, k);
		for (int j = 1; j <= pcnt && 1ll * i * prime[j] <= maxn; ++ j) {
			is_prime[i * prime[j]] = false; if (!(i % prime[j])) break;
		}
	}
}

int id1[N], id2[N]; ui val[N * 2], ptot[N * 2], d, all;

#define id(x) (x <= d ? id1[x] : id2[all / (x)])

void Min25_Sieve(ui n) {
	int cnt = 0;
	for (ui i = 1; i <= n; i = n / (n / i) + 1)
		val[id(n / i) = ++ cnt] = n / i, ptot[cnt] = val[cnt] - 1;

	for (int i = 1; i <= pcnt && 1ll * prime[i] * prime[i] <= n; ++ i)
		for (int j = 1; j <= cnt && 1ll * prime[i] * prime[i] <= val[j]; ++ j)
			ptot[j] -= ptot[id(val[j] / prime[i])] - (i - 1);
}

ui S(ui n, int cur) {
	if (n <= 1 || (ui)prime[cur] > n) return 0;
	ui res = (ptot[id(n)] - (cur - 1)) * Pow[cur - 1];
	for (int i = cur; i <= pcnt && 1ll * prime[i] * prime[i] <= n; ++ i) {
		ui prod = prime[i];
		for (int e = 1; 1ll * prod * prime[i] <= n; ++ e, prod *= prime[i])
			res += S(n / prod, i + 1) + Pow[i];
	}
	return res;
}

bitset<N * 2> vis; ui M[N * 2];

ui Calc(ui n) {
	if (vis[id(n)]) return M[id(n)];
	ui res = S(n, 1) + ptot[id(n)];
	for (ui i = 2, ni; i <= n; i = ni + 1)
		ni = n / (n / i), res -= (ni - i + 1) * Calc(n / i);
	return vis[id(n)] = true, M[id(n)] = res;
}

int main () {

	File();

	ui n = read(); k = read();

	Linear_Sieve(d = sqrt(n + .5));

	all = n; Min25_Sieve(n); 

	ui ans = 0, Last = 0;
	for (ui i = 1, ni; i <= n; i = ni + 1) {
		ni = n / (n / i); ui res = Calc(ni);
		ans += (n / i) * (n / i) * (res - Last); Last = res;
	}
	printf ("%u\n", ans);

	return 0;

}

__EOF__

本文作者：zjp_shadow
本文链接：https://www.cnblogs.com/zjp-shadow/p/10202492.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！