用Dirichlet卷积计算互质幂和 fₖ(n)

hbghlyj

给出$n=\prod_{i=1}^w p_i^{a_i}$，$p_i$都为质数且互不相同，$w,a_i\ge1$。求
$$
f_k(n)=\sum_{\substack{1\leqslant i\leqslant n\\\gcd(i,n)=1}}i^k
$$
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}$显然是积性函数，所以对每个质因子求出之后乘起来即可。提前高斯消元出来所有$a_i$。

1. #include <algorithm>
2. #include <cstdio>
3. #include <cstring>
4. const int D = 105;
5. const int mod = 1000000007;
6. typedef long long LL;
7. inline LL pow_mod(LL a, LL b) {
8.   LL ans = 1;
9.   for (a %= mod, (b += mod - 1) %= mod - 1; b; b >>= 1, a = a * a % mod)
10.     if (b & 1) ans = ans * a % mod;
11.   return ans;
12. }
13. inline LL inv(LL a) { return pow_mod(a, -1); }
14. int d;
15. LL a[D];
16. void solve() {
17.   static LL A[D][D];
18.   LL t = 0;
19.   for (int i = 0; i <= d; ++i) {
20.     LL j = 1;
21.     for (int k = 0; k <= d; ++k) A[i][k] = j = j * (i + 1) % mod;
22.     A[i][d + 1] = ((t += d ? A[i][d - 1] : 1) %= mod);
23.   }
24.   for (int i = 0; i <= d; ++i) {
25.     int j = i;
26.     while (!A[j][i]) ++j;
27.     for (int k = i; k <= d + 1; ++k) std::swap(A[i][k], A[j][k]);
28.     LL inv1 = inv(A[i][i]);
29.     for (int j = i; j <= d + 1; ++j)
30.       A[i][j] = A[i][j] * inv1 % mod;
31.     for (int j = i + 1; j <= d; ++j)
32.       for (int k = d + 1; k >= i; --k)
33.         A[j][k] = (A[j][k] - A[j][i] * A[i][k] % mod) % mod;
34.   }
35.   for (int i = d; ~i; --i) {
36.     a[i + 1] = A[i][d + 1];
37.     for (int j = i - 1; ~j; --j)
38.       A[j][d + 1] = (A[j][d + 1] - A[j][i] * A[i][d + 1] % mod) % mod;
39.   }
40. }
41. const int N = 1050;
42. int p[N], q[N];
43. int main() {
44.   int w;
45.   scanf("%d%d", &d, &w);
46.   solve();
47.   LL ans = 0, n = 1;
48.   for (int i = 0; i < w; ++i) {
49.     scanf("%d%d", &p[i], &q[i]);
50.     n = n * pow_mod(p[i], q[i]) % mod;
51.   }
52.   LL y = 1;
53.   for (int i = 1; i <= d + 1; ++i) {
54.     y = y * n % mod;
55.     LL t = a[i] * y % mod;
56.     for (int j = 0; j < w; ++j)
57.       t = t * (1 - pow_mod(p[j], d - i)) % mod;
58.     ans = (ans + t) % mod;
59.   }
60.   printf("%d\n", (int)((ans + mod) % mod));
61.   return 0;
62. }

Copy the Code