『题解』Luogu-P6271 [湖北省队互测2014]一个人的数论
Description
-
给出非负整数 \(d\) 和正整数 \(n\) 的质因数分解式 \(n = \prod_{i = 1}^{\omega} p_i^{\alpha_i}\)，请求出
\[\left(f_d(n) = \sum_{i = 1}^n [\gcd(i, n) = 1]\, i^d \right) \bmod (10^9 + 7) \] -
\(0\le d\le 100, 1\le \omega\le 10^3, 2\le p_i\le 10^9, 1\le \alpha_i \le 10^9\)。
Solution
\[\begin{aligned}
f_d(n)
& = \sum_{i = 1}^n [\gcd(i, n) = 1] \, i^d \\
& = \sum_{i = 1}^n i^d \sum_{k\mid \gcd(i, n)} \mu(k) \\
& = \sum_{k \mid n} \mu(k) \sum_{i = 1}^n [k\mid i] \, i^d \\
& = \sum_{k \mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k}} i^d
\end{aligned}
\]
后面的等幂和直接用伯努利数做就可以了。
我们知道
\[\begin{aligned}
S_m(n)
& = \sum_{k = 0}^{n - 1} k^m \\
& = \dfrac{1}{m + 1} \sum_{k = 0}^m \dbinom{m + 1}{k} B_k n^{m + 1 - k}
\end{aligned}
\]
设出系数
\[f_i = \dfrac{1}{d + 1} \dbinom{d + 1}{i} B_i
\]
系数可以 \(\Omicron(d^2)\) 预处理。
则
\[\begin{aligned}
f_d(n)
& = \sum_{k\mid n} \mu(k) k^d S_d\left(\dfrac{n}{k} + 1 \right) \\
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 0}^d f_i \left(\dfrac{n}{k} + 1 \right)^{d + 1 - i}
\end{aligned}
\]
你发现这个 \(\left(\dfrac{n}{k} + 1 \right)^{m + 1 - i}\) 没法处理，但如果是 \(\left(\dfrac{n}{k}\right)^{m + 1 - i}\) 就很好了——它可以和前面的 \(k^d\) 相乘。
退回到上一步
\[\begin{aligned}
f_d(n)
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k}} i^d \\
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d + \sum_{k\mid n} \mu(k) k^d \left(\dfrac{n}{k}\right)^d \\
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d + n^d \sum_{k\mid n} \mu(k) \\
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d + n^d [n = 1] \\
\end{aligned}
\]
观察数据范围
\(1\le \omega\le 10^3, 2\le p_i\le 10^9, 1\le \alpha_i \le 10^9\)。
说明 \(n\ne 1\)。
则
\[\begin{aligned}
f_d(n)
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 1}^{\frac{n}{k} - 1} i^d \\
& = \sum_{k \mid n} \mu(k) k^d S_d\left(\dfrac{n}{k} \right) \\
& = \sum_{k\mid n} \mu(k) k^d \sum_{i = 0}^d f_i \left(\dfrac{n}{k}\right)^{d + 1 - i} \\
& = \sum_{k\mid n} \mu(k) \sum_{i = 0}^d f_i n^{d + 1 - i} k^{i - 1} \\
& = \sum_{i = 0}^d f_i n^{d + 1 - i} \sum_{k\mid n} \mu(k) k^{i - 1}
\end{aligned}
\]
你会发现
\[\sum_{k\mid n} \mu(k) k^{i - 1} = [(\mu\cdot \operatorname{Id}_{i - 1}) * \mathbf{1}](n)
\]
所以这是一个积性函数，但是也不可能线性筛啊。
我们设
\[\sum_{k\mid n} \mu(k) k^{i - 1} = g_i(n)
\]
因为 \(g_i(n)\) 是积性函数，所以可以按质因数拆开：
\[g_i(n) = \prod_{k = 1}^{\omega} g_i(p_k^{\alpha_k})
\]
则只用考虑质因数幂的情况。
你会发现
\[g_i(p^{\alpha}) = \sum_{k\mid p^{\alpha}} \mu(k) k^{i - 1}
\]
其中只有 \(k = 1\) 和 \(k = p\) 时 \(\mu(k)\) 才有贡献。
而
\[\begin{cases}
\mu(1) 1^{i - 1} = 1 \\
\mu(p) p^{i - 1} = - p^{i - 1}
\end{cases}
\]
所以
\[g_i(p^{\alpha}) = 1 - p^{i - 1}
\]
再贴一遍式子
\[f_d(n) = \sum_{i = 0}^d f_i n^{d + 1 - i} \sum_{k\mid n} \mu(k) k^{i - 1}
\]
计算
\[g_i(n) = \prod_{k = 1}^{\omega} 1 - p_k^{i - 1}
\]
是 \(\Omicron(\omega \log i)\) 的，求个和就是 \(\Omicron(\omega d \log d)\)，但其实可以把 \(p^{i - 1}\) 往上递推做到 \(\Omicron(d\omega)\)。
预处理伯努利数 \(f_i\) 是 \(\Omicron(d^2)\) 的，而 \(n^{d + 1 - i}\) 是 \(\Omicron(d\log d)\) 的，远不及前面的 \(\Omicron(d\omega)\)。
综上，时间复杂度为 \(\Omicron(d^2 + d\omega)\)。
Code
注意预处理要处理到 \(d + 1\)。
//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;
const int MOD = 1e9 + 7;
const int MAXD = 105;
int add(int a, int b) {return (a + b) % MOD;}
int sub(int a, int b) {return (a - b + MOD) % MOD;}
int mul(int a, int b) {return (ll)a * b % MOD;}
int qpow(int a, int b) {int base = a, ans = 1; while (b) {if (b & 1) ans = mul(ans, base); base = mul(base, base); b >>= 1;} return ans;}
int GetInv(int a) {return qpow(a, MOD - 2);}
struct Bernoulli
{
int fac[MAXD], fac_inv[MAXD], inv[MAXD];
void pre(int d)
{
fac[0] = fac_inv[0] = inv[1] = 1;
for (int i = 1; i <= d; i++)
{
fac[i] = mul(fac[i - 1], i);
}
fac_inv[d] = GetInv(fac[d]);
for (int i = d - 1; i >= 1; i--)
{
fac_inv[i] = mul(fac_inv[i + 1], i + 1);
inv[i + 1] = mul(fac_inv[i + 1], fac[i]);
}
}
int C(int n, int m)
{
return mul(mul(fac[n], fac_inv[n - m]), fac_inv[m]);
}
int B[MAXD], f[MAXD];
void cal(int d)
{
B[0] = 1;
f[0] = inv[d + 1];
for (int m = 1; m <= d; m++)
{
int sum = 0;
for (int k = 0; k < m; k++)
{
sum = add(sum, mul(C(m + 1, k), B[k]));
}
B[m] = mul(MOD - sum, inv[m + 1]);
f[m] = mul(mul(inv[d + 1], C(d + 1, m)), B[m]);
}
}
}Ber;
const int MAXW = 1005;
int p[MAXW], a[MAXW], xsy[MAXW];
int main()
{
int d, w;
scanf("%d%d", &d, &w);
Ber.pre(d + 1), Ber.cal(d);
int n = 1;
for (int i = 1; i <= w; i++)
{
scanf("%d%d", p + i, a + i);
n = mul(n, qpow(p[i], a[i]));
xsy[i] = GetInv(p[i]);
}
int ans = 0;
for (int i = 0; i <= d; i++)
{
int pro = 1;
for (int j = 1; j <= w; j++)
{
pro = mul(pro, sub(1, xsy[j]));
xsy[j] = mul(xsy[j], p[j]);
}
ans = add(ans, mul(mul(Ber.f[i], qpow(n, d + 1 - i)), pro));
}
printf("%d\n", ans);
return 0;
}
