「学习笔记」单位根反演
单位根反演
公式推导
\[[ n \mid k ] = \frac{1}{n}\sum\limits_{ i = 0 }^{ n - 1 } \omega_{ n }^{ ik }
\]
证明：
当 \(n\mid k\) 时，\(\omega_{n}^{ik}=1\)，显然成立；
当 \(n\nmid k\) 时，对右边项进行等比数列求和得到 \(\frac{1}{n}\frac{\omega_{n}^{nk}-1}{\omega_{n}^{k}-1}=0\)。
单位根反演常用于求某个多项式特定倍数的系数和。
\[\begin{aligned}
\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}[x^{ik}]f(x)&=\sum_{i=0}^{n}[k\mid i][x^i]f(x)\\
&=\sum_{i=0}^{n}[x^i]f(x)\frac{1}{k}\sum_{j=0}^{k-1}\omega_{k}^{ji}\\
&=\frac{1}{k}\sum_{i=0}^{n}a_i\sum_{j=0}^{k-1}\omega_{k}^{ij}\\
&=\frac{1}{k}\sum_{i=0}^{n}\sum_{j=0}^{k-1}a_i(\omega_{k}^{j})^i\\
&=\frac{1}{k}\sum_{j=0}^{k-1}f(\omega_{k}^{j})
\end{aligned}
\]
同理可推出
\[\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}[x^{ik+r}]f(x)=\frac{1}{k}\sum_{i=0}^{k-1}f(\omega_{k}^{i})\omega_{k}^{-ir}\\
\]
例题
\(\color{blue}{\text{[LOJ6485] LJJ } 学二项式定理}\)
求 \(\big{[}\sum\limits_{i=0}^{n}\big{(}\binom{n}{i}\cdot s^{i}\cdot a_{i\operatorname{mod} 4}\big{)} \big{]}\operatorname{mod} 998244353\)。
\(1\le T\le10^5,1\le n\le 10^{18},1\le s,a_0,a_1,a_2,a_3\le 10^8\)。
\[\begin{aligned}
\sum\limits_{i=0}^{n}\big{(}\binom{n}{i}\cdot s^{i}\cdot a_{i\operatorname{mod} 4}\big{)}&=\sum_{j=0}^{3}\sum_{i=0}^{n}\binom{n}{i}\cdot s^{i}\cdot a_j[4|i-j]\\
&=\sum_{j=0}^{3}\sum_{i=0}^{n}\binom{n}{i}\cdot s^i
\cdot a_j\cdot\frac{1}{4}\sum_{k=0}^{3}\omega_{4}^{k(i-j)}\\
&=\frac{1}{4}\sum_{k=0}^{3}\sum_{j=0}^{3}a_j\cdot\omega_{4}^{-jk}\sum_{i=0}^{n}\binom{n}{i}\cdot s^i\cdot \omega_{4}^{ik}\\
&=\frac{1}{4}\sum_{k=0}^{3}\sum_{j=0}^{3}a_j\cdot\omega_{4}^{-jk}(s\omega_{4}^{k}+1)^n
\end{aligned}
\]
直接计算即可，时间复杂度 \(O(\log n)\)。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 998244353;
int T, s, ans, a[5]; ll n;
inline ll qpow(ll a, ll b) {
ll ans = 1;
for (; b; b >>= 1, a = a * a % mod) if (b & 1) ans = ans * a % mod;
return ans;
}
const int omega = qpow(3, (mod - 1) / 4);
const int omega_ = qpow(omega, mod - 2);
const int inv4 = qpow(4, mod - 2);
int main() {
scanf("%d", &T);
while (T --) {
int ans = 0;
scanf("%lld%d", &n, &s);
for (int i = 0; i < 4; ++ i) scanf("%d", a + i);
for (int k = 0; k < 4; ++ k)
for (int j = 0; j < 4; ++ j)
(ans += a[j] * qpow(omega_, j * k) % mod * (qpow((s * qpow(omega, k) + 1) % mod, n)) % mod) %= mod;
printf("%d\n", (ll)ans * inv4 % mod);
}
return 0;
}
\([\color{blue}{\text{AtCoder Regular Contest 145 F - Modulo Sum of Increasing Sequences]}}\)
对于每个 \(K=0,1,\ldots,\text{MOD}-1\)，计数长度为 \(N\)，值域为 \([0,M]\)，总和对 \(\text{MOD}\) 取模后结果为 \(K\) 的不降序列 \(A\) 的数量，答案对 \(998244353\) 取模。
\(1\le N,M\le10^6,1\le\text{MOD}\le 500\)。
令 \(b_i=a_i+i-1\)，则有 \(\forall i<n, b_i<b_{i+1}\)。问题转化成了在 \([0,M+N-1]\) 中选取 \(N\) 个数使得总和对 \(\text{MOD}\) 取模后结果为 \(K\)。下令 \(n=M+N\)。为了方便，记 \(p=\text{MOD}\)。
先不考虑长度的限制，列出 OGF
\[F(x)=\prod_{i=0}^{n-1}(1+x^i)
\]
答案即为
\[\frac{1}{p}\sum_{i=0}^{p-1}F(\omega_{p}^{i})\omega_{p}^{-ik}=\frac{1}{p}\sum_{i=0}^{p-1}\prod_{j=0}^{n-1}(1+\omega_{p}^{ij})\omega_{p}^{-ik}
\]
\(\color{blue}{\textbf{Special Case: } p\mid n}\)
假设 \(p\mid n\)，记 \(d=\gcd (p,i)\)，则
\[\frac{1}{p}\sum_{i=0}^{p-1}\prod_{j=0}^{n-1}(1+\omega_{p}^{ij})\omega_{p}^{-ik}=\frac{1}{p}\sum_{i=0}^{p-1}(\prod_{j=0}^{\frac{p}{d}-1}(1+\omega_{\frac{p}{d}}^{ij})^{d})^{\frac{n}{p}}\omega_{p}^{-ik}
\]
用分圆多项式作简化
\[x^n-1=\prod_{i=0}^{n-1}(x-\omega_{n}^{i})
\]
带入 \(x=-1\)，得
\[(-1)^{n}-1=(-1)^n\prod_{i=0}^{n-1}(1+\omega_{n}^{i})
\]
所以
\[\prod_{i=0}^{n-1}(1+\omega_{n}^{i})=1-(-1)^{n}
\]
那么
\[\begin{aligned}
\frac{1}{p}\sum_{i=0}^{p-1}(\prod_{j=0}^{\frac{p}{d}-1}(1+\omega_{\frac{p}{d}}^{ij})^{d})^{\frac{n}{p}}\omega_{p}^{-ik}&=\frac{1}{p}\sum_{i=0}^{p-1}\big{(}\prod_{j=0}^{\frac{p}{d}-1}(1+\omega_{\frac{p}{d}}^{ij}) \big{)}^\frac{dn}{p}\omega_{p}^{-ik}\\
&=\frac{1}{p}\sum_{i=0}^{p-1}(1-(-1)^{\frac{p}{d}} )^\frac{dn}{p}\omega_{p}^{-ik}\\
&=\frac{1}{p}\sum_{d\mid p}(1-(-1)^{\frac{p}{d}})^\frac{dn}{p}(\sum_{\gcd(i,p)=d}\omega_{p}^{-ik})\\
&=\frac{1}{p}\sum_{d\mid p}(1-(-1)^{\frac{p}{d}})^\frac{dn}{p}(\sum_{\gcd(i,\frac{p}{d})=1}\omega_{\frac{p}{d}}^{-ik})\\
&=\frac{1}{p}\sum_{d\mid p}(1-(-1)^{\frac{p}{d}})^\frac{dn}{p}(\sum_{\gcd(i,\frac{p}{d})=1}\omega_{\frac{p}{d}}^{-ik})\\
&=\frac{1}{p}\sum_{d\mid p}(1-(-1)^{\frac{p}{d}})^\frac{dn}{p}(\sum_{j\mid\frac{p}{d}}\mu(j)\sum_{i=0}^{\frac{p}{dj}-1}\omega_{\frac{p}{dj}}^{-ik})\\
&=\frac{1}{p}\sum_{d\mid p}(1-(-1)^{\frac{p}{d}})^\frac{dn}{p}(\sum_{j\mid\frac{p}{d}}\mu(j)[\frac{p}{dj}\mid k]\frac{p}{dj})
\end{aligned}
\]
现在加上长度的限制，在初始的 OGF 中多加一个元 \(y\)，变为
\[\prod_{i=0}^{n-1}(1+x^iy)
\]
推导后的结果为
\[[y^{N}]\frac{1}{p}\sum_{d\mid p}(1-(-y)^{\frac{p}{d}})^\frac{dn}{p}(\sum_{j\mid\frac{p}{d}}\mu(j)[\frac{p}{dj}\mid k]\frac{p}{dj})
\]
可以用二项式定理 \(O(p^2)\) 求出。
\(\color{blue}{\textbf{General Case: } p\nmid n}\)
根据上面的推导，可以 \(O(p^2)\) 求出 \(p\mid n\) 的情况，我们在此基础上做 \(p\nmid n\) 的情况。
首先令 \(n'=\lfloor\frac{n}{p}\rfloor\cdot p\)，\(O(p^2)\) 求出 \(0\sim n'-1\) 的答案。
考虑添加进 \(n'\sim n-1\) 之间的数，也就是求出在 \(n'\sim n-1\) 之间选取 \(i\) 个数满足这些数的和 \(\operatorname{mod}{p}=j\) 的方案数，最后再和之前的答案进行合并。
因为 \(n'-n<p\)，所以这部分直接用 \(dp\) 解决，时间复杂度为 \(O(p^3)\)。
至此，本题解决。总时间复杂度 \(O(N+M+p^3)\)。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 505, M = 2e6 + 5, mod = 998244353;
int n, m, p, _n, mu[N], ans[N], dp[N][N];
ll fac[M], inv[M], Inv[M];
vector<int> divisors[N];
vector<pair<int, vector<int>>> gr;
inline void precalc(int n) {
fac[0] = inv[0] = Inv[0] = fac[1] = inv[1] = Inv[1] = 1;
for (int i = 2; i < n; ++ i) {
fac[i] = fac[i - 1] * i % mod;
Inv[i] = (mod - mod / i) * Inv[mod % i] % mod;
inv[i] = inv[i - 1] * Inv[i] % mod;
}
}
inline ll C(int n, int m) {
if (n < m || n < 0 || m < 0) return 0;
return fac[n] * inv[m] % mod * inv[n - m] % mod;
}
int main() {
precalc(M);
scanf("%d%d%d", &n, &m, &p), n += m, m = n - m, _n = n / p * p;
dp[0][0] = 1;
for (int i = 1; i <= n - _n; ++ i)
for (int k = i; k; -- k)
for (int j = 0; j < p; ++ j)
(dp[k][j] += dp[k - 1][(j - (i + _n - 1) % p + p) % p]) %= mod;
for (int i = 1; i <= p; ++ i)
for (int j = i; j <= p; j += i)
divisors[j].emplace_back(i);
mu[1] = 1;
for (int i = 1; i <= p; ++ i)
for (int j = i + i; j <= p; j += i)
mu[j] -= mu[i];
for (int i = 2; i <= p; ++ i) (mu[i] += mod) %= mod;
for (auto d : divisors[p]) {
vector<int> c(p);
for (int k = 0; k < p; ++ k)
for (auto j : divisors[p / d])
if ((k * d * j) % p == 0)
(c[k] += (ll)p / d / j * mu[j] % mod) %= mod;
gr.emplace_back(d, c);
}
for (int use = 0; use <= min(n - _n, m); ++ use) {
int leave = m - use;
vector<int> coef(p);
for (auto &&[d, c] : gr) {
int b = p / d; ll co;
if (leave % b > 0) continue;
if (b & 1) co = Inv[p];
else co = (leave / b) & 1 ? mod - Inv[p] : Inv[p];
(co *= C(d * _n / p, leave / b)) %= mod;
for (int i = 0; i < p; ++ i)
(coef[i] += co * c[i] % mod) %= mod;
}
for (int i = 0; i < p; ++ i)
for (int j = 0; j < p; ++ j)
(ans[(i + j) % p] += (ll)coef[i] * dp[use][j] % mod) %= mod;
}
int move = (m - 1LL) * m / 2 % p;
for (int i = 0; i < p; ++ i) printf("%d ", ans[(i + move) % p]);
return 0;
}
