Chirp Z-transform
Similar to the Discrete Fourier Transform, the Chirp Z-transform is an algorithm that, given a polynomial \(f(x) = \sum_{i = 0}^{m - 1} f_i x^i \in \mathbb{C}\lbrack x\rbrack\) and \(q \in \mathbb{C} \setminus \{0\}\) , computes \(f(1), f(q), \dots, f(q^{n - 1})\) , without requiring \(q\) to be a root of unity. It can also be used for number-theoretic transforms. The following sections will introduce the Chirp Z-transform and its inverse.
By definition, the Chirp Z-transform can be written as
\[
\operatorname{\mathsf{CZT}}_n : \left(f(x), q\right) \mapsto
\begin{bmatrix}
f(1) & f(q) & \cdots & f\left(q^{n - 1}\right)
\end{bmatrix}
\]
where \(f(x) := \sum_{i = 0}^{m - 1} f_i x^i \in \mathbb{C}\lbrack x\rbrack\) and \(q \in \mathbb{C} \setminus \{0\}\) .
Bluestein's Algorithm
Consider
\[
ij = \binom{i}{2} + \binom{-j}{2} - \binom{i - j}{2}
\]
where \(i, j \in \mathbb{Z}\) . We can construct
\[
\begin{aligned}
G(x) & := \sum_{i = -(m - 1)}^{n - 1} q^{-\binom{i}{2}}x^i, \\
F(x) & := \sum_{i = 0}^{m - 1} f_i q^{\binom{-i}{2}}x^i.
\end{aligned}
\]
Here \(G(x) \in \mathbb{C}\left\lbrack x, x^{-1}\right\rbrack\) , and for \(i = 0, \dots, n - 1\) we have
\[
\begin{aligned}
\left\lbrack x^i\right\rbrack\left(G(x)F(x)\right) &=
\sum_{j = 0}^{m - 1}\left(\left(\left\lbrack x^{i - j}\right\rbrack G(x)\right)\left(\left\lbrack x^j\right\rbrack F(x)\right)\right) \\
&= \sum_{j = 0}^{m - 1} f_j q^{\binom{-j}{2} - \binom{i - j}{2}} \\
&= q^{-\binom{i}{2}} f\left(q^i\right)
\end{aligned}
\]
and \(q^{\binom{i + 1}{2}} = q^{\binom{i}{2}}\cdot q^i\) , \(\binom{-i}{2} = \binom{i + 1}{2}\) . This can be computed via a single polynomial multiplication. This algorithm is known as Bluestein's algorithm.
Template (P6800【Template】Chirp Z-Transform )
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 std :: vector < uint > CZT ( const std :: vector < uint >& f , uint q , int n ) {
if ( f . empty () || n == 0 ) return std :: vector < uint > ( n );
const int m = f . size ();
if ( q == 0 ) {
std :: vector < uint > res ( n , f [ 0 ]);
for ( int i = 1 ; i < m ; ++ i )
if (( res [ 0 ] += f [ i ]) >= MOD ) res [ 0 ] -= MOD ;
return res ;
}
// H[i] = q^(binom(i + 1, 2))
std :: vector < uint > H ( std :: max ( m , n - 1 ));
H [ 0 ] = 1 ; // H[0] = q^0
uint qi = 1 ; // qi = q^i
for ( int i = 1 ; i < ( int ) H . size (); ++ i ) {
qi = ( ull ) qi * q % MOD ;
// H[i] = q^(binom(i, 2)) * q^i
H [ i ] = ( ull ) H [ i - 1 ] * qi % MOD ;
}
// F[i] = f[i] * q^(binom(i + 1, 2))
std :: vector < uint > F ( m );
for ( int i = 0 ; i < m ; ++ i ) F [ i ] = ( ull ) f [ i ] * H [ i ] % MOD ;
std :: vector < uint > G_p ( m + n - 1 );
// G[i] = q^(-binom(i, 2)), -(m - 1) ≤ i < n
uint * const G = G_p . data () + ( m - 1 );
const uint iq = InvMod ( q );
// G[-(m - 1)] = q^(-binom(-(m - 1), 2)),
// binom(-(m - 1), 2) = (-(m - 1)) * (-(m - 1) - 1) / 2
// = (m - 1) * m / 2
G [ - ( m - 1 )] = PowMod ( iq , ( ull )( m - 1 ) * m / 2 );
uint qmi = PowMod ( q , m - 1 ); // q^(-i), -(m - 1) ≤ i < n
for ( int i = - ( m - 1 ) + 1 ; i < n ; ++ i ) {
// q^(-binom(i, 2)) = q^(-binom(i - 1, 2)) * q^(-(i - 1))
G [ i ] = ( ull ) G [ i - 1 ] * qmi % MOD ;
// q^(-i) = q^(-(i - 1)) * q^(-1)
qmi = ( ull ) qmi * iq % MOD ;
}
// res[i] = q^(-binom(i, 2)) * f(q^i), 0 ≤ i < n
std :: vector < uint > res = MiddleProduct ( G_p , F );
for ( int i = 1 ; i < n ; ++ i ) res [ i ] = ( ull ) res [ i ] * H [ i - 1 ] % MOD ;
return res ;
}
The inverse Chirp Z-transform can be written as
\[
\operatorname{\mathsf{ICZT}}_n :
\left(
\begin{bmatrix} f(1) & f(q) & \cdots & f\left(q^{n - 1}\right)
\end{bmatrix},q
\right)
\mapsto f(x)
\]
where \(f(x) \in \mathbb{C}\left\lbrack x\right\rbrack_{< n}\) and \(q \in \mathbb{C} \setminus \{0\}\) , and \(q^i \neq q^j\) holds for all \(i \neq j\) , which is the condition for polynomial interpolation.
Bostan–Schost Algorithm
Recall the Lagrange interpolation formula is
\[
f(x) = \sum_{i = 0}^{n - 1}\left(f\left(x_i\right)\prod_{0 \leq j < n \atop j \neq i} \frac{x - x_j}{x_i - x_j}\right)
\]
and \(x_i \neq x_j\) for all \(i \neq j\) . As in Fast Polynomial Interpolation , let \(M(x) := \prod_{i = 0}^{n - 1}\left(x - x_i\right)\) . By L'Hôpital's rule,
\[
M'(x_i) = \lim_{x \to x_i} \frac{M(x)}{x - x_i} = \prod_{0 \leq j < n \atop j \neq i}\left(x_i - x_j\right)
\]
The modified Lagrange interpolation formula is
\[
f(x) = M(x)\left(\sum_{i = 0}^{n - 1}\frac{f\left(x_i\right)/M'(x_i)}{x - x_i}\right)
\]
Thus we now have
\[
f(x) = M(x)\left(\sum_{i = 0}^{n - 1}\frac{f\left(q^i\right)/M'\left(q^i\right)}{x - q^i}\right)
\]
where \(M(x)=\prod_{j = 0}^{n - 1}\left(x - q^j\right)\) . If we set \(n\) to be even, let \(n = 2k\) and \(H(x) := \prod_{j = 0}^{k - 1}\left(x - q^j\right)\) , then
\[
M(x) = H(x) \cdot q^{k^2} \cdot H\left(\frac{x}{q^k}\right)
\]
This allows us to compute \(M(x)\) quickly. Then use Bluestein's algorithm to compute \(M'(1), \dots, M'(q^{n - 1})\) . Let \(c_i := f\left(q^i\right)/M'\left(q^i\right)\) . We have
\[
f(x) = M(x)\left(\sum_{i = 0}^{n - 1}\frac{c_i}{x - q^i}\right)
\]
Since \(\deg f(x) < n\) , we only need to compute \(\sum_{i = 0}^{n - 1}\frac{c_i}{x - q^i}\bmod{x^n}\) , where \(\frac{c_i}{x - q^i} \in \mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack\) , i.e.,
\[
\begin{aligned}
\sum_{i = 0}^{n - 1}\frac{c_i}{x - q^i} \bmod x^n &=
-\sum_{i = 0}^{n - 1}\left(\sum_{j = 0}^{n - 1} c_i q^{-i(j+1)}x^j\right) \\
&= -\sum_{j = 0}^{n - 1} C\left(q^{-j - 1}\right) x^j
\end{aligned}
\]
where \(C(x) = \sum_{i = 0}^{n - 1} c_i x^i\) . We can use Bluestein's algorithm to compute \(C\left(q^{-1}\right), \dots, C\left(q^{-n}\right)\) .
In summary, we perform the following computations:
Compute \(M(x)\) using the divide-and-conquer method;
Compute \(M'(1), \dots, M'(q^{n - 1})\) using Bluestein's algorithm;
Compute \(C\left(q^{-1}\right), \dots, C\left(q^{-n}\right)\) using Bluestein's algorithm;
Compute \(f(x)\) via a single polynomial multiplication.
The time complexity of each step is equal to that of multiplying two polynomials of degree at most \(n\) .
Template Implementation
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 std :: vector < uint > InvCZT ( const std :: vector < uint >& f , uint q ) {
if ( f . empty ()) return {};
const int n = f . size ();
if ( q == 0 ) {
// f[0] = f(1), f[1] = f(0)
assert ( n <= 2 );
if ( n == 1 ) return { f [ 0 ]}; // deg(f) < 1
return { f [ 1 ], ( f [ 0 ] + MOD - f [ 1 ]) % MOD }; // deg(f) < 2
}
// prod[0 ≤ i < n] (x - q^i)
const auto DaC = [ q ]( auto && DaC , int n ) -> std :: vector < uint > {
if ( n == 1 ) return { MOD - 1 , 1u };
// H = prod[0 ≤ i < ⌊n/2⌋] (x - q^i)
const std :: vector < uint > H = DaC ( DaC , n / 2 );
// HH = H(x / q^(⌊n/2⌋)) * q^(⌊n/2⌋^2)
std :: vector < uint > HH = H ;
const uint iqn = InvMod ( PowMod ( q , n / 2 ));
uint qq = PowMod ( q , ( ull )( n / 2 ) * ( n / 2 ));
for ( int i = 0 ; i <= n / 2 ; ++ i )
HH [ i ] = ( ull ) HH [ i ] * qq % MOD , qq = ( ull ) qq * iqn % MOD ;
std :: vector < uint > res = Product ( H , HH );
if ( n & 1 ) {
res . resize ( n + 1 );
const uint qnm1 = MOD - PowMod ( q , n - 1 );
for ( int i = n ; i > 0 ; -- i ) {
if (( res [ i ] += res [ i - 1 ]) >= MOD ) res [ i ] -= MOD ;
res [ i - 1 ] = ( ull ) res [ i - 1 ] * qnm1 % MOD ;
}
}
return res ;
};
const std :: vector < uint > M = DaC ( DaC , n );
// C[i] = (M'(q^i))^(-1)
std :: vector < uint > C = InvModBatch ( CZT ( Deriv ( M ), q , n ));
// C[i] = f(q^i) / M'(q^i)
for ( int i = 0 ; i < n ; ++ i ) C [ i ] = ( ull ) f [ i ] * C [ i ] % MOD ;
// C(x) ← C(q^(-1)x)
const uint iq = InvMod ( q );
uint qmi = 1 ;
for ( int i = 0 ; i < n ; ++ i )
C [ i ] = ( ull ) C [ i ] * qmi % MOD , qmi = ( ull ) qmi * iq % MOD ;
C = CZT ( C , iq , n );
for ( int i = 0 ; i < n ; ++ i )
if ( C [ i ] != 0 ) C [ i ] = MOD - C [ i ];
std :: vector < uint > res = Product ( M , C );
res . resize ( n );
return res ;
}
References
Bostan, A. (2010). Fast algorithms for polynomials and matrices. JNCF 2010. Algorithms Project, INRIA.
更新历史 在 GitHub 上编辑此页！ OI-wiki CC BY-SA 4.0 和 SATA