AIOJ

Like 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})$ ; it does not require $q$ to be a root of unity. It can also be used for number-theoretic transforms. The Chirp Z transform and its inverse are introduced below.

Chirp Z transform

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 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}

Also, $q^{\binom{i + 1}{2}} = q^{\binom{i}{2}}\cdot q^i$ and $\binom{-i}{2} = \binom{i + 1}{2}$ . The computation can be completed with one polynomial multiplication; this algorithm is called the Bluestein algorithm.

Template (P6800 [Template] Chirp Z-Transform)

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;
}

Inverse Chirp Z transform

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$ for all $i \neq j$ . This is the condition for polynomial interpolation.

Bostan–Schost algorithm

Recall that 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)

where $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)

Therefore,

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 $n$ is 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)

which lets us compute $M(x)$ quickly. We then use the Bluestein algorithm to compute $M'(1), \dots, M'(q^{n - 1})$ . Let $c_i := f\left(q^i\right)/M'\left(q^i\right)$ . Then

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$ . That is,

\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 the Bluestein algorithm to compute $C\left(q^{-1}\right), \dots, C\left(q^{-n}\right)$ .

In short, we perform the following computations:

Compute $M(x)$ using decrease and conquer.
Use the Bluestein algorithm to compute $M'(1), \dots, M'(q^{n - 1})$ .
Use the Bluestein algorithm to compute $C\left(q^{-1}\right), \dots, C\left(q^{-n}\right)$ .
Compute $f(x)$ with one polynomial multiplication.

The time complexity of every step equals the time complexity of multiplying two polynomials of degree at most $n$ .

Template implementation

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.