【学习笔记】Chirp-Z Transform

其实很简单（）

前置知识：FFT / NTT

001. P6800 【模板】Chirp Z-Transform

给定一个 \(n\) 项多项式 \(P(x)\) 以及 \(c, m\)，请计算 \(P(c^0),P(c^1),\dots,P(c^{m-1})\)。所有答案都对 \(998244353\) 取模。

题解：

\(P(c^i)\) 的值为 \(\sum\limits_{j=0}^{n-1}a_jc^{ij}\)，因此要求的东西其实就是：

\(\forall i\in[0,m),F_i=\sum\limits_{j=0}^{n-1}a_jc^{ij}\bmod 998244353\)

这是一个二维卷积的形式，考虑使用 Chirp-Z Transform 优化。有核心公式：

\[ij=\binom{i+j}2-\binom i2-\binom j2 \]

用这个东西来替换上式的 \(ij\)，可得：

\[\begin{aligned} F_i&=\sum\limits_{j=0}^{n-1}a_jc^{ij}\\&=\sum\limits_{j=0}^{n-1}a^jc^{\binom{i+j}2-\binom i2-\binom j2}\\&=c^{-\binom i2}\sum\limits_{j=0}^{n-1}a^jc^{\binom{i+j}2-\binom j2}\\&={\large{c^{-\binom i2}}}\sum\limits_{j=0}^{n-1}\left(\left(a^jc^{-\binom j2}\right)\left(c^\binom{i+j}2\right)\right) \end{aligned} \]

设 \(A_i=a_ic^{-\binom i2},B_i=c^{\binom i2}\)，那么上述式子形如一个等和卷积，翻转 \(B\) 得到 \(B'\)，设 \(C=A\odot B\) 即 \(A\) 和 \(B\) 的卷积，即可求得后半部分的和式。而前面的 \(c^{-\binom i2}\) 可以单 \(\log\) 计算也可以线性递推。总时间复杂度为 \(O(n\log n)\)，瓶颈在于 NTT。

我的实现细节处理的不怎么好，但是也能过（）

namespace Luminescent
{struct DSU{int fa[N];inline DSU() { iota(fa, fa + N, 0); }inline int find(int x) { return x == fa[x] ? x : fa[x] = find(fa[x]); }inline int merge(int x, int y){x = find(x), y = find(y);if (x != y)return fa[x] = y, 1;return 0;}};inline void add(int &x, int v) { x = (x + v) % mod; }inline void sub(int &x, int v) { x = (x - v + mod) % mod; }inline int power(int a, int b, int c){int ans = 1;while (b){if (b & 1)ans = ans * a % c;a = a * a % c, b >>= 1;}return ans;}inline int inversion(int x) { return power(x, mod - 2, mod); }inline int varphi(int x){int phi = 1;for (int i = 2; i * i <= x; ++i)if (x % i == 0){phi *= (i - 1);x /= i;while (x % i == 0)phi *= i, x /= i;}if (x > 1)phi *= (x - 1);return phi;}
} using namespace Luminescent;namespace Poly
{const int g = 3;int rev[N * 8];void ntt(int *a, int n, int mode){int bit = 1;while ((1 << bit) < n)++bit;for (int i = 0; i < n; ++i){rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));if (i < rev[i])swap(a[i], a[rev[i]]);}for (int l = 2; l <= n; l <<= 1){int x = power(g, (mod - 1) / l, mod);if (mode == 1)x = inversion(x);for (int i = 0; i < n; i += l){int v = 1;for (int j = 0; j < l / 2; ++j, v = v * x % mod){int v1 = a[i + j], v2 = a[i + j + l / 2] * v % mod;a[i + j] = (v1 + v2) % mod, a[i + j + l / 2] = (v1 - v2 + mod) % mod;}}}}// calc convolution: c[i] = \sum\limits_{j=0}^i (a[j] * b[i - j])void convolution(int *a, int n, int *b, int m, int *c){int tn = n, tm = m;n = n + m + 2;while (__builtin_popcount(n) > 1)++n;// cerr << "n = " << n << '\n';for (int i = tn + 1; i <= n + 1; ++i)a[i] = 0;for (int i = tm + 1; i <= n + 1; ++i)b[i] = 0;ntt(a, n, 0), ntt(b, n, 0);for (int i = 0; i < n; ++i)c[i] = a[i] * b[i] % mod;ntt(c, n, 1);const int inv_n = inversion(n);for (int i = 0; i <= n + m; ++i)c[i] = c[i] * inv_n % mod;}
}namespace Loyalty
{int a[N], A[N * 8], B[N * 8], C[N * 8];inline void main([[maybe_unused]] int _ca, [[maybe_unused]] int atc){int n, c, m;cin >> n >> c >> m;for (int i = 0; i < n; ++i)cin >> a[i];for (int i = 0; i < max(n, m); ++i)A[i] = a[i] * inversion(power(c, i * (i - 1) / 2, mod)) % mod;for (int i = 0; i < max(n, m) + max(n, m); ++i)B[max(n, m) + max(n, m) - i - 1] = power(c, i * (i - 1) / 2, mod);// cerr << "???\n";Poly::convolution(A, max(n, m), B, max(n, m) + max(n, m), C);// cerr << "!!!\n";for (int j = 0; j < m; ++j){int sum = 0;sum = C[max(n, m) + max(n, m) - j - 1];// for (int i = 0; i < n; ++i)//     sum = (sum + A[i] * B[n + n - i - j - 1] % mod) % mod;cout << sum * inversion(power(c, j * (j - 1) / 2, mod)) % mod << ' ';}cout << '\n';}
}signed main()
{// freopen("b.in", "r", stdin);// freopen("b.oo", "w", stdout);cin.tie(0)->sync_with_stdio(false);cout << fixed << setprecision(15);int T = 1;// cin >> T;for (int ca = 1; ca <= T; ++ca)Loyalty::main(ca, T);return 0;
}

002. P14561 [CXOI2025] 我常常追忆过去

对于整数 \(n\)，有一张 \(n\) 个点的竞赛图，图上每一条边 \((u,v)\)（\(u<v\)）的方向有 \(\frac12\) 的概率是由 \(u\) 连向 \(v\) 的，另外 \(\frac12\) 的概率是由 \(v\) 连向 \(u\) 的。

现在给定 \(m\)，对于每个 \(n=1\ldots m\)，你需要求出 \(n\) 个点的竞赛图上强连通分量的数目的期望，答案对 \(998244353\) 取模。

这么会夹带私货

首先先思考，假设给定了 \(n\) 的值，怎么快速计算答案。

容易证明若 \(A,B\) 两个点集之间，不存在任何一条边 \(u\to v\) 使得 \(u\in B,v\in A\)，那么 \(A,B\) 之间就形成了一个新的 SCC。

考虑直接计数这样被分割的 SCC 数量。考虑枚举集合 \(A\) 的大小 \(i\)，那么有 \(\binom ni\) 种选择 \(A\) 的方法。定义 \(B\) 集合是 \(A\) 在全集下的补集，则若不存在 \(B\) 连向 \(A\) 的边，则两个点集之间的全部 \(2^{|A|\times|B|}\) 条边都必须由 \(A\) 集合连向 \(B\) 集合即已被定向，因此答案即为 \(1+\sum\limits_{i=1}^{n-1}\binom ni2^{-i(n-i)}\)，简单细节处理或使用 \(O(\sqrt n)-O(1)\) 的光速幂均可做到 \(O(n)\) 计算答案。

然后考虑优化到快速计算 \(1\sim n\) 的答案。记竞赛图大小为 \(n\) 时的期望为 \(E_n\)，则有：

\[E_n=1+\sum\limits_{i=1}^{n-1}\binom ni2^{-i(n-i)}=\sum\limits_{i=1}^n\binom ni2^{-i(n-i)} \]

给这个东西搞一个 EGF 然后只保留系数，可以得到：

\[\frac{E_n}{n!}=\sum\limits_{i=1}^n\frac{2^{-i(n-i)}}{i!(n-i)!} \]

记 \(p\equiv 2^{-1}\pmod{998244353}\)，则上式可以改写为：

\[\frac{E_n}{n!}=\sum\limits_{i=1}^n\frac{p^{i(n-i)}}{i!(n-i)!} \]

然后特殊处理掉边界上 \(i=0\) 的部分：

\[\frac{E_n}{n!}+\frac1{n!}=\sum\limits_{i=0}^n\frac{p^{i(n-i)}}{i!(n-i)!} \]

这是一个二维卷积，可以用 Chirp-Z 变换和 NTT 将其优化至 \(O(n\log n)\) 求解。

具体而言：考虑经典套路，有 \(i(n-i)=\binom n2-\binom i2-\binom{n-i}2\)，于是开始推柿子：

\[\begin{aligned}&\sum\limits_{i=0}^n\frac{p^{i(n-i)}}{i!(n-i)!}\\ =&\sum\limits_{i=0}^n\frac{p^{\binom n2-\binom i2-\binom{n-i}2}}{i!(n-i)!}\\ =&\sum\limits_{i=0}^n\frac{p^{\binom n2}\times p^{-\binom i2}\times p^{-\binom{n-i}2}}{i!(n-i)!}\\ =&p^{\binom n2}\sum\limits_{i=0}^n\left(\frac{p^{-\binom i2}}{i!}\times\frac{p^{-\binom{n-i}2}}{(n-i)!}\right)\\ \end{aligned} \]

设 \(A_i=\dfrac{p^{-\binom i2}}{i!}\)，则设 \(C=A\odot A\) 即自己和自己的卷积，上述表达式的值可以记作 \(p^{\binom n2}C_n\)。因为 \(998244353\) 是 NTT 友好模数，所以可以 \(O(n\log n)\) 求解上述问题。

namespace Luminescent
{struct DSU{int fa[N];inline DSU() { iota(fa, fa + N, 0); }inline int find(int x) { return x == fa[x] ? x : fa[x] = find(fa[x]); }inline int merge(int x, int y){x = find(x), y = find(y);if (x != y)return fa[x] = y, 1;return 0;}};inline void add(int &x, int v) { x = (x + v) % mod; }inline void sub(int &x, int v) { x = (x - v + mod) % mod; }inline int power(int a, long long b, int c){int sum = 1;while (b){if (b & 1)sum = 1ll * sum * a % c;a = 1ll * a * a % c, b >>= 1;}return sum;}inline int inversion(int x) { return power(x, mod - 2, mod); }inline int varphi(int x){int phi = 1;for (int i = 2; i * i <= x; ++i)if (x % i == 0){phi *= (i - 1);x /= i;while (x % i == 0)phi *= i, x /= i;}if (x > 1)phi *= (x - 1);return phi;}
} using namespace Luminescent;namespace Poly
{const int g = 3;int rev[N * 4];void ntt(int *a, int n, int mode){int bit = 1;while ((1 << bit) < n)++bit;for (int i = 0; i < n; ++i){rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));if (i < rev[i])swap(a[i], a[rev[i]]);}for (int l = 2; l <= n; l <<= 1){int x = power(g, (mod - 1) / l, mod);if (mode == 1)x = inversion(x);for (int i = 0; i < n; i += l){int v = 1;for (int j = 0; j < l / 2; ++j, v = v * x % mod){int v1 = a[i + j], v2 = a[i + j + l / 2] * v % mod;a[i + j] = (v1 + v2) % mod, a[i + j + l / 2] = (v1 - v2 + mod) % mod;}}}}// calc convolution: c[i] = \sum\limits_{j=0}^i (a[j] * b[i - j])void convolution(int *a, int n, int *b, int m, int *c){int tn = n, tm = m;n = n + m + 2;while (__builtin_popcount(n) > 1)++n;// cerr << "n = " << n << '\n';for (int i = tn + 1; i <= n + 1; ++i)a[i] = 0;for (int i = tm + 1; i <= n + 1; ++i)b[i] = 0;ntt(a, n, 0), ntt(b, n, 0);for (int i = 0; i < n; ++i)c[i] = a[i] * b[i] % mod;ntt(c, n, 1);const int inv_n = inversion(n);for (int i = 0; i <= n + m; ++i)c[i] = c[i] * inv_n % mod;}
}namespace Loyalty
{int fac[N], inv[N], ifac[N];int A[N << 2], B[N << 2], C[N << 2];inline void main([[maybe_unused]] int _ca, [[maybe_unused]] int atc){int m;cin >> m;for (int i = 0; i < 2; ++i)fac[i] = inv[i] = ifac[i] = 1;for (int i = 2; i < N; ++i){fac[i] = fac[i - 1] * i % mod;inv[i] = mod - inv[mod % i] * (mod / i) % mod;ifac[i] = ifac[i - 1] * inv[i] % mod;}const int inv_2 = mod + 1 >> 1;for (int i = 0; i <= m; ++i)A[i] = B[i] = inversion(power(inv_2, i * (i - 1) / 2, mod)) * ifac[i] % mod;Poly::convolution(A, m, B, m, C);for (int i = 1; i <= m; ++i)cout << (C[i] * power(inv_2, i * (i - 1) / 2, mod) % mod * fac[i] - 1 + mod) % mod << '\n';}
}

003. P5293 [HNOI2019] 白兔之舞

前置知识：矩阵优化 dp，单位根反演

先考虑 \(L\) 很小怎么做。容易想到枚举一共走了 \(i\) 步，后面的方案数形如组合数乘以矩阵的形式。设转移矩阵为 \(M\)，则对于每个 \(t\in[0,k)\cap\mathbb N\)，答案可以表示为：

\[\begin{aligned}&\sum\limits_{i\equiv t\pmod k}\binom Li(M_i)_{x,y}\\ =&\sum\limits_i\binom Li(M_i)_{x,y}[i\equiv t\bmod k]\\ =&\sum\limits_i\binom Li(M_i)_{x,y}[k\mid i-t]\\ =&\sum\limits_i[k\mid i-t]\binom Li(M_i)_{x,y}\\ =&\sum\limits_i\frac1k\sum\limits_{j=0}^{k-1}\omega_k^{j(i-t)}\binom Li(M_i)_{x,y}\\ =&\frac1k\sum\limits_{i}\sum\limits_{j=0}^{k-1}\omega_k^{j(i-t)}\binom Li(M_i)_{x,y}\\ =&\frac1k\sum\limits_i\sum\limits_{j=0}^{k-1}\omega_k^{ij}\times\omega_k^{-tj}\binom Li(M_i)_{x,y}\\ =&\frac1k\sum\limits_{j=0}^{k-1}\omega_k^{-tj}\sum\limits_{i}\omega_k^{ij}\binom Li(M_i)_{x,y}\\ =&\frac1k\sum\limits_{j=0}^{k-1}\omega_k^{-tj}\sum\limits_i\binom Li\textrm{I}^{L-i}\omega_k^{ij}(M_i)_{x,y}\\ =&\frac 1k\sum\limits_{j=0}^{k-1}\omega_k^{-tj}((\textrm{I}+\omega_k^j\times M)^L)_{x,y} \end{aligned} \]

把上述式子写成关于 \(\omega_k^j\) 的函数，即设 \(f(x)=\frac 1k\sum\limits_{j=0}^{k-1}x^j((\mathrm{I}+\omega_k^j\times M)^L)_{x,y}\)，问题即为对每个 \(t\in[0,k)\cap\mathbb N\) 求 \(f(\omega_k^{-t})\) 的值。注意到 \(f(x)\) 是一个多项式函数，而要求的值形如二维卷积的形式，因此可以使用 Chirp-Z Transform 优化。

因为这个题模数不固定，所以需要使用 MTT，这里使用了拆系数的 FFT 来计算答案（拆系数是因为卡精度）。因为这个题太卡精度了，所以这里使用了 long double。

namespace Luminescent
{const double pi = acos(-1);const ld pi_l = acosl(-1);struct DSU{int fa[N];inline DSU() { iota(fa, fa + N, 0); }inline int find(int x) { return x == fa[x] ? x : fa[x] = find(fa[x]); }inline int merge(int x, int y){x = find(x), y = find(y);if (x != y)return fa[x] = y, 1;return 0;}};inline void add(int &x, int v) { x = (x + v) % mod; }inline void sub(int &x, int v) { x = (x - v + mod) % mod; }inline int power(int a, int b, int c){int sum = 1;while (b){if (b & 1)sum = 1ll * sum * a % c;a = 1ll * a * a % c, b >>= 1;}return sum;}inline int inversion(int x, int mod) { return power(x, mod - 2, mod); }inline int varphi(int x){int phi = 1;for (int i = 2; i * i <= x; ++i)if (x % i == 0){phi *= (i - 1);x /= i;while (x % i == 0)phi *= i, x /= i;}if (x > 1)phi *= (x - 1);return phi;}
} using namespace Luminescent;int n;namespace LA_Matrix
{struct Matrix{int a[3][3];inline Matrix() { memset(a, 0, sizeof a); }inline Matrix(int val){for (int i = 0; i < n; ++i)for (int j = 0; j < n; ++j)a[i][j] = val;}};const inline Matrix I(){Matrix o;o.a[0][0] = o.a[1][1] = o.a[2][2] = 1;return o;}inline Matrix operator+(const Matrix &l, const Matrix &r){Matrix res;for (int i = 0; i < n; ++i)for (int j = 0; j < n; ++j)res.a[i][j] = (l.a[i][j] + r.a[i][j]) % mod;return res;}inline Matrix operator*(const Matrix &l, const Matrix &r){Matrix res;for (int k = 0; k < n; ++k)for (int i = 0; i < n; ++i)for (int j = 0; j < n; ++j)res.a[i][j] = (res.a[i][j] + l.a[i][k] * r.a[k][j] % mod) % mod;return res;}inline Matrix operator*(const Matrix &l, int r){r %= mod;Matrix res;for (int i = 0; i < n; ++i)for (int j = 0; j < n; ++j)res.a[i][j] = l.a[i][j] * r % mod;return res;}inline Matrix power(Matrix a, int b){if (!b)return I();Matrix res = power(a, b >> 1);res = res * res;if (b & 1)res = res * a;return res;}
}namespace Poly
{int rev[N];using cauto = complex<ld>;inline void fft(cauto *a, int n, int mode){int bit = 1;while ((1 << bit) < n)++bit;for (int i = 0; i < n; ++i)if (i < rev[i])swap(a[i], a[rev[i]]);for (int len = 2; len <= n; len <<= 1){cauto x = cauto(cosl(pi_l * 2. / len), sinl(pi_l * 2. / len));if (mode == 1)x = cauto(cosl(pi_l * 2. / len), -sinl(pi_l * 2. / len));for (int i = 0; i < n; i += len){cauto v = cauto(1, 0);for (int j = 0; j < len / 2; ++j, v = v * x){cauto v1 = a[i + j], v2 = a[i + j + len / 2] * v;a[i + j] = v1 + v2, a[i + j + len / 2] = v1 - v2;}}}}cauto A1[N], B1[N], A2[N], B2[N], C1[N], C2[N], C3[N];inline void convolution(int *a, int *b, int len, int *c){int n = len << 1;while (n & (n - 1))++n;int bit = 1;while ((1 << bit) < n)++bit;for (int i = 0; i < n; ++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));for (int i = 0; i < n; ++i)A1[i].real(a[i] >> 15), A2[i].real(a[i] & 32767);for (int i = 0; i < n; ++i)B1[i].real(b[i] >> 15), B2[i].real(b[i] & 32767);fft(A1, n, 1), fft(A2, n, 1), fft(B1, n, 1), fft(B2, n, 1);for (int i = 0; i < n; ++i)C1[i] = A2[i] * B2[i], C2[i] = A1[i] * B2[i] + A2[i] * B1[i], C3[i] = A1[i] * B1[i];fft(C1, n, -1), fft(C2, n, -1), fft(C3, n, -1);for (int i = 0; i < n; ++i)C1[i].real(C1[i].real() / n), C2[i].real(C2[i].real() / n), C3[i].real(C3[i].real() / n);for (int i = 0; i < n; ++i){int v0 = (int)(C1[i].real() + 0.5) % mod;int v1 = (int)(C2[i].real() + 0.5) % mod;int v2 = (int)(C3[i].real() + 0.5) % mod;c[i] = (v2 * (1ll << 30) % mod + v1 * (1ll << 15) % mod + v0) % mod;}}
}namespace Loyalty
{vector<int> fj;inline int chk(int x, int p){for (int &i : fj)if (power(x, (p - 1) / i, p) == 1)return 0;return 1;}inline int find_g(int p){for (int i = 2; ; ++i)if (chk(i, p))return i;return -1;}int k, L, x, y;int m1[N], m2[N], m3[N];inline void main([[maybe_unused]] int _ca, [[maybe_unused]] int atc){cin >> n >> k >> L >> x >> y >> mod;int _ = mod - 1;for (int i = 2; i * i <= _; ++i)if (_ % i == 0){fj.emplace_back(i);while (_ % i == 0)_ /= i;}if (_ > 1)fj.emplace_back(_);const int g = find_g(mod), omega = power(g, (mod - 1) / k, mod), inv_k = inversion(k, mod);// cerr << omega << '\n';LA_Matrix::Matrix M_map;for (int i = 0; i < n; ++i)for (int j = 0; j < n; ++j)cin >> M_map.a[i][j];// cerr << "m1: ";for (int i = 0; i < k; ++i)m1[i] = (LA_Matrix::power(LA_Matrix::I() + M_map * power(omega, i, mod), L)).a[x - 1][y - 1] * power(omega, i * (i - 1) / 2, mod) % mod;//, cerr << m1[i] << ' ';// cerr << '\n';// return;for (int i = 0; i < k + k; ++i)m2[k + k - i - 1] = power(omega, mod - i * (i - 1) / 2 % (mod - 1) - 1, mod);//, cerr << m2[k + k - i - 1] << ' ';// return;// cerr << '\n';// for (int i = 0; i < k; ++i)//     for (int j = 0; j < k + k; ++j)//         m3[i + j] = (m3[i + j] + m1[i] * m2[j] % mod) % mod;Poly::convolution(m1, m2, k + k, m3);// cerr << "m3: ";// for (int i = 0; i < k + k; ++i)//     cerr << m3[i] << ' ';// cerr << '\n';for (int i = 0; i < k; ++i)cout << m3[k + k - i - 1] * power(omega, i * (i - 1) / 2, mod) % mod * inv_k % mod << '\n';}
}signed main()
{// freopen("4.in", "r", stdin);// freopen("test.out1", "w", stdout);cin.tie(0)->sync_with_stdio(false);cout << fixed << setprecision(15);int T = 1;// cin >> T;for (int ca = 1; ca <= T; ++ca)Loyalty::main(ca, T);return 0;
}

004. CF1054H Epic Convolution

咕咕咕，不知道什么时候能更新（）