Loj 572 Misaka Network 与求和

Posted on | Views:

杜教筛 + 魔改 min_25 筛.

为了方便,记 $f(i)$ 就表示次大质因子的 $k​$ 次幂.

$$
\begin{aligned}
ans&=\sum_{i=1}^n\sum_{j=1}^n f(\gcd(i,j)) \\
&=\sum_{d=1}^n f(d) \sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{^{\lfloor\frac n d\rfloor}} [\gcd(i,j)=1] \\
&=\sum_{d=1}^n \lfloor\frac{n}{d}\rfloor^2\sum_{t|d}\mu(t)f(\frac{d}{t})\\
&=\sum_{d=1}^n \lfloor\frac{n}{d}\rfloor^2\cdot (f*\mu)(d)
\end{aligned}
$$

对 $d$ 整除分块,就需要多次询问 $f* \mu$ 的前缀和.

记 $h(n)=\sum_{i=1}^n (f* \mu)(i)$ .

用杜教筛的套路,可以得出

$$
\sum_{i=1}^n f(n) = \sum_{i=1}^n (f \times \mu \times I) (i) = \sum_{i=1}^n h(\lfloor \frac n i \rfloor) \\
h(n)=\sum_{i=1}^n f(i)-\sum_{i=2}^n h(\lfloor \frac n i \rfloor)
$$

Latex 渲染有点问题,就用 $\times$ 代替 $*$ 作为狄利克雷卷积的符号了.

计算 $h$ 时整除分块,递归下去计算.

每次还会询问 $f$ 的前缀和,可以把 min_25 筛魔改一下,在提出因子时,如果它是次大因子,才计入贡献.

此外,还要把 $\le n$ 的每个质数贡献的 $1$ 加入 $h(n)$ .

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
//%std
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int uint;
inline uint read()
{
	uint out = 0, fh = 1;
	char jp = getchar();
	while ((jp > '9' || jp < '0') && jp != '-')
		jp = getchar();
	if (jp == '-')
		fh = -1, jp = getchar();
	while (jp >= '0' && jp <= '9')
		out = out * 10 + jp - '0', jp = getchar();
	return out * fh;
}
uint fpow(uint a, int b)
{
	uint res = 1;
	while (b)
	{
		if (b & 1)
			res = res * a;
		a = a * a;
		b >>= 1;
	}
	return res;
}
const int N = 2e5 + 10;
int n, k, sqr, w[N], tot = 0, id1[N], id2[N];
int cnt = 0, ism[N], prime[N]; 
uint pk[N], g[N], _h[N];
int vis_h[N];
int id(int x)
{
	if (x <= sqr)
		return id1[x];
	return id2[n / x];
}
void init()
{
	ism[1] = 1, pk[1] = 1;
	for (int i = 2; i <= sqr; ++i)
	{
		if (!ism[i])
			prime[++cnt] = i, pk[i] = fpow(i, k);
		for (int j = 1; j <= cnt && i * prime[j] <= sqr; ++j)
		{
			ism[i * prime[j]] = 1;
			if (i % prime[j] == 0)
				break;
			pk[i * prime[j]] = pk[i] * pk[prime[j]];
		}
	}
	for (int l = 1, r; l <= n; l = r + 1)
	{
		r = n / (n / l);
		w[++tot] = n / l;
		if (n / l <= sqr)
			id1[n / l] = tot;
		else
			id2[n / (n / l)] = tot;
	}
	for (int i = 1; i <= tot; ++i)
		g[i] = w[i] - 1;
	for (int j = 1; j <= cnt; ++j)
		for (int i = 1; i <= tot && prime[j] * prime[j] <= w[i]; ++i)
		{
			int k = id(w[i] / prime[j]);
			g[i] -= g[k] - j + 1;
		}
}
uint S(int x, int j)
{
	if (x <= 1 || prime[j] > x)
		return 0;
	uint res = 0;
	for (int i = j; i <= cnt && prime[i] * prime[i] <= x; ++i)
		for (ll pw = prime[i]; pw * prime[i] <= x; pw *= prime[i])
			res += S(x / pw, i + 1) + (g[id(x / pw)] - i + 1) * pk[prime[i]];
	return res;
}
uint h(int x)
{
	int pos = id(x);
	if (vis_h[pos])
		return _h[pos];
	vis_h[pos] = 1;
	uint res = S(x, 1) + g[pos];
	for (int l = 2, r; l <= x; l = r + 1)
	{
		r = x / (x / l);
		res -= h(x / l) * (uint)(r - l + 1);
	}
	return _h[pos] = res; 
}
int main()
{
	n = read(), k = read();
	sqr = ceil(sqrt(n));
	init();
	uint ans = 0, lst = 0, cur;
	for (int l = 1, r; l <= n; l = r + 1)
	{
		r = n / (n / l);
		cur = h(r);
		ans += (cur - lst) * (uint)(n / l) * (uint)(n / l);
		lst = cur;
	}
	printf("%u\n", ans);
	return 0;
}