奇怪的数学题 | 51nod1847

题目链接

给定 \(n,k\),求 \(\sum\limits_{i=1}^n\sum\limits_{j=1}^n\text{sgcd}(i,j)^k\)

其中 \(\text{sgcd}(i,j)\) 表示 \(i,j\) 的次大公约数。特别地,如果 \(\text{gcd}(i,j)=1\),则 \(\text{sgcd}(i,j)=0\)

答案对 \(2^{32}\) 取模。

\(n \le 10^9,k \le 50\)

\(p_i\) 表示第 \(i\) 个质数,\(d_i\) 表示 \(i\) 的最小质因子。

\(\text{sgcd}(i,j)=\dfrac{\gcd(i,j)}{d_{\gcd(i,j)}}\)

考虑枚举 \(\gcd\)\[ \begin{aligned} &\sum_{i=1}^n\sum_{j=1}^n\text{sgcd}(i,j)^k\\ &=\sum_{c=2}^n(\frac c{d_c})^k\sum_{i=1}^{\lfloor\frac nc\rfloor}\sum_{j=1}^{\lfloor\frac nc\rfloor}[gcd(i,j)=1]\\ &=\sum_{c=2}^n(\frac c{d_c})^k(2\sum_{i=1}^{\lfloor\frac nc\rfloor}\varphi(i)-1) \end{aligned} \] 现在的问题是算 \((\frac x{d_x})^k\)\(\varphi(x)\)\(\lfloor \frac nc \rfloor\) 处的前缀和,后者直接杜教筛即可。

对于前者,由于涉及到 \(d_x\) 考虑 Min-25 筛

\[ g(n,i) = \sum_{j=2}^n [j \in P \lor d_j > p_i]j^k\\ f(n,i) = \sum_{j=2}^n [j \not\in P \land d_j \le p_i](\frac j{d_j})^k\\ h(h,i) = \sum_{j=2}^n [j \in P \lor d_j > p_i] \] 有递推 \[ g(n, i) = g(n, i - 1) - p_i^k(g(\lfloor\frac n {p_i} \rfloor, i - 1) - g(p_i-1,i))\\ f(n,i) = f(n,i-1)+g(\lfloor\frac n {p_i} \rfloor,i-1) - g(p_i-1,i)\\ h(n,i)=h(n,i-1)-h(\lfloor\frac n {p_i} \rfloor, i - 1) + h(p_i-1,i) \] 初始化 \[ g(n,0)=\sum_{i=2}^ni^k\\ f(n,0)=0\\ h(n,0)=n-1 \] 这里需要求自然数等幂和

\(n < p_i^2\) 时都不需要转移,因此这是一个标准的 Min-25 筛。

复杂度 \(O(\frac {n^{\frac 34}}{\log n}+n^{\frac 23})\)

代码:

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
#include <bits/stdc++.h>
#define rep(i, l, r) for(int i = (l); i <= (r); i++)
#define per(i, r, l) for(int i = (r); i >= (l); i--)
#define mem(a, b) memset(a, b, sizeof a)
#define For(i, l, r) for(int i = (l), i##e = (r); i < i##e; i++)
#define pb push_back

using namespace std;
typedef unsigned int U;

int n, K;
namespace Sum {
U S[55][55];
void pre() {
S[0][0] = 1;
rep(i, 1, K) rep(j, 1, i) S[i][j] = S[i - 1][j - 1] + (U)j * S[i - 1][j];
}
U qry(int n, U re = 0) {
rep(i, 1, K) {
U t = 1;
rep(j, n + 1 - i, n + 1) if(j % (i + 1)) t *= j;
else t *= j / (i + 1);
re += t * S[K][i];
}
return re - 1u;
}
};
namespace Du {
int m, pid, prm[100000];
U phi[1000005], S[1005];
void pre() {
m = pow(n, 2. / 3);
phi[1] = 1;
rep(i, 2, m) {
if(!phi[i]) phi[i] = i - 1, prm[++pid] = i;
for(int j = 1; i * prm[j] <= m; j++)
if(i % prm[j]) phi[i * prm[j]] = phi[i] * phi[prm[j]];
else { phi[i * prm[j]] = phi[i] * prm[j]; break; }
}
rep(i, 1, m) phi[i] += phi[i - 1];
}
U qry(int i) {
if(i <= m) return phi[i];
if(S[n / i]) return S[n / i];
U res = i * (i + 1ll) / 2;
for(int l = 2, r; l <= i; l = r + 1) {
r = i / (i / l);
res -= qry(i / l) * (r - l + 1);
}
return S[n / i] = res;
}
};
namespace M25 {
constexpr int N = sqrt(1e9) + 5;
int m;
U g1[N], g2[N], f1[N], f2[N], h1[N], h2[N];
void pre() {
m = sqrt(n);
rep(i, 1, m) {
g1[i] = Sum::qry(i), g2[i] = Sum::qry(n / i);
h1[i] = i - 1, h2[i] = n / i - 1;
}
rep(p, 2, m) if(h1[p] ^ h1[p - 1]) {
int w1 = m / p, w3 = p * p, w2 = min(m, n / w3);
int j, d = n / p;
U x = 1, gx = g1[p - 1], hx = h1[p - 1];
rep(i, 1, K) x *= p;
rep(i, 1, w1) {
j = i * p, h2[i] -= h2[j] - hx;
f2[i] += g2[j] - gx, g2[i] -= x * (g2[j] - gx);
}
rep(i, w1 + 1, w2) {
j = d / i, h2[i] -= h1[j] - hx;
f2[i] += g1[j] - gx, g2[i] -= x * (g1[j] - gx);
}
per(i, m, w3) {
j = i / p, h1[i] -= h1[j] - hx;
f1[i] += g1[j] - gx, g1[i] -= x * (g1[j] - gx);
}
}
}
U qry(int i) {return i <= m ? f1[i] + h1[i] : f2[n / i] + h2[n / i]; }
};
int main() {
cin >> n >> K;
Sum::pre(), Du::pre(), M25::pre();
U as = 0;
for(int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
as += (2u * Du::qry(n / l) - 1u) * (M25::qry(r) - M25::qry(l - 1));
}
cout << as;
return 0;
}