BZOJ3601 一个人的数论

Description

定义 \[ f_k(n)=\sum_{\substack{1\leqslant i\leqslant n\\gcd(i,n)=1}}i^k \]

给出\(n=\prod_{i=1}^w p_i^{a_i}\),求\(f_k(n)\)\(1\leqslant w\leqslant 1000, 1\leqslant q_i,a_i\leqslant 10^9\)。保证\(p_i\)都为质数且互不相同。

Solution

\(g_k(n)=\sum_{i=1}^ni^k\),则 \[ \begin{aligned} g_k(n)&=\sum_{i=1}^ni^k\\ &=\sum_{d|n}\sum_{\substack{1\leqslant i\leqslant n\\gcd(i,n)=d}}i^k\\ &=\sum_{d|n}d^k\sum_{\substack{1\leqslant i\leqslant \frac nd\\gcd\left(i,\frac nd\right)=d}}i^k\\ &=\sum_{d|n}d^kf_k\left(\frac nd\right) \end{aligned} \] 也即\(g_k(n) = (n^k) * f_k(n)\)\(*\)表示狄利克雷卷积)。

由于\((n^k)*(\mu(n)n^k)=[n=1]\),所以\(f_k(n)=(\mu(n)n^k)*g_k(n)\)

显然\(g_k(n)\)可以写成\(\sum_{i=1}^{k+1}a_in^i\),那么 \[ \begin{aligned} f_k(n)&=\sum_{d|n}\mu(d)d^k\sum_{i=1}^{k+1}a_i\left(\frac nd\right)^i\\ &=\sum_{i=1}^{k+1}a_i\sum_{d|n}\mu(d)d^k\left(\frac nd\right)^i\\ &=\sum_{i=1}^{k+1}a_in^i\sum_{d|n}\mu(d)d^{k-i}\\ \end{aligned} \] \(\sum_{d|n}\mu(d)d^{k-i}\)显然是积性函数,所以对每个质因子\(O(1)\)求出之后乘起来即可。

所有\(a_i\)提前高消出来就行了。

Code

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
#include <algorithm>
#include <cstdio>
#include <cstring>
const int D = 105;
const int mod = 1000000007;
typedef long long LL;
inline LL pow_mod(LL a, LL b) {
LL ans = 1;
for (a %= mod, (b += mod - 1) %= mod - 1; b; b >>= 1, a = a * a % mod)
if (b & 1) ans = ans * a % mod;
return ans;
}
inline LL inv(LL a) { return pow_mod(a, -1); }
int d;
LL a[D];
void solve() {
static LL A[D][D];
LL t = 0;
for (int i = 0; i <= d; ++i) {
LL j = 1;
for (int k = 0; k <= d; ++k) A[i][k] = j = j * (i + 1) % mod;
A[i][d + 1] = ((t += d ? A[i][d - 1] : 1) %= mod);
}
for (int i = 0; i <= d; ++i) {
int j = i;
while (!A[j][i]) ++j;
for (int k = i; k <= d + 1; ++k) std::swap(A[i][k], A[j][k]);
LL inv1 = inv(A[i][i]);
for (int j = i; j <= d + 1; ++j)
A[i][j] = A[i][j] * inv1 % mod;
for (int j = i + 1; j <= d; ++j)
for (int k = d + 1; k >= i; --k)
A[j][k] = (A[j][k] - A[j][i] * A[i][k] % mod) % mod;
}
for (int i = d; ~i; --i) {
a[i + 1] = A[i][d + 1];
for (int j = i - 1; ~j; --j)
A[j][d + 1] = (A[j][d + 1] - A[j][i] * A[i][d + 1] % mod) % mod;
}
}
const int N = 1050;
int p[N], q[N];
int main() {
int w;
scanf("%d%d", &d, &w);
solve();
LL ans = 0, n = 1;
for (int i = 0; i < w; ++i) {
scanf("%d%d", &p[i], &q[i]);
n = n * pow_mod(p[i], q[i]) % mod;
}
LL y = 1;
for (int i = 1; i <= d + 1; ++i) {
y = y * n % mod;
LL t = a[i] * y % mod;
for (int j = 0; j < w; ++j)
t = t * (1 - pow_mod(p[j], d - i)) % mod;
ans = (ans + t) % mod;
}
printf("%d\n", (int)((ans + mod) % mod));
return 0;
}