Chirp Z 变换

与离散傅里叶变换类似，Chirp Z 变换是给出多项式 $f(x)={\sum }_{i=0}^{m-1}{f}_i{x}^i\in \mathbb{C}[x]$ 和 $q\in \mathbb{C}\setminus \{0\}$ 求出 $f(1),f(q),\dots ,f(q^{n-1})$ 的一种算法，不要求 $q$ 为单位根。也可用于数论变换。后文将介绍 Chirp Z 变换与其逆变换。

Chirp Z 变换

根据定义，Chirp Z 变换可以写作

$${\mathrm{CZT}}_n:\left(f(x),q\right)\mapsto \left[\begin{array}{cccc}f(1)& f(q)& \cdots & f\left(q^{n-1}\right)\end{array}\right]$$

其中 $f(x):={\sum }_{i=0}^{m-1}{f}_i{x}^i\in \mathbb{C}[x]$ 且 $q\in \mathbb{C}\setminus \{0\}$。

Bluestein 算法

考虑

$$ij=\left(\genfrac{}{}{0ex}{}{i}{2}\right)+\left(\genfrac{}{}{0ex}{}{-j}{2}\right)-\left(\genfrac{}{}{0ex}{}{i-j}{2}\right)$$

其中 $i,j\in \mathbb{Z}$，我们可以构造

$$\begin{array}{rl}G(x)& :=\sum _{i=-(m-1)}^{n-1}{q}^{-\left(\genfrac{}{}{0ex}{}{i}{2}\right)}{x}^i,\\ F(x)& :=\sum _{i=0}^{m-1}{f}_i{q}^{\left(\genfrac{}{}{0ex}{}{-i}{2}\right)}{x}^i.\end{array}$$

其中 $G(x)\in \mathbb{C}\left[x,x^{-1}\right]$，且对于 $i=0,\dots ,n-1$ 我们有

$$\begin{array}{rl}\left[x^i\right]\left(G(x)F(x)\right)& =\sum _{j=0}^{m-1}\left(\left(\left[x^{i-j}\right]G(x)\right)\left(\left[x^j\right]F(x)\right)\right)\\ & =\sum _{j=0}^{m-1}{f}_j{q}^{\left(\genfrac{}{}{0ex}{}{-j}{2}\right)-\left(\genfrac{}{}{0ex}{}{i-j}{2}\right)}\\ & =q^{-\left(\genfrac{}{}{0ex}{}{i}{2}\right)}f\left(q^i\right)\end{array}$$

且 $q^{\left(\genfrac{}{}{0ex}{}{i+1}{2}\right)}=q^{\left(\genfrac{}{}{0ex}{}{i}{2}\right)}\cdot q^i$，$\left(\genfrac{}{}{0ex}{}{-i}{2}\right)=\left(\genfrac{}{}{0ex}{}{i+1}{2}\right)$。可以由一次多项式乘法完成求算，该算法被称为 Bluestein 算法。

模板（ P6800【模板】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;
}

逆 Chirp Z 变换

逆 Chirp Z 变换可以写作

$${\mathrm{ICZT}}_n:\left(\left[\begin{array}{cccc}f(1)& f(q)& \cdots & f\left(q^{n-1}\right)\end{array}\right],q\right)\mapsto f(x)$$

其中 $f(x)\in \mathbb{C}{\left[x\right]}_{<n}$ 且 $q\in \mathbb{C}\setminus \{0\}$，并且 $q^i\ne q^j$ 对于所有 $i\ne j$ 成立，这是多项式插值的条件。

Bostan–Schost 算法

回顾 Lagrange 插值公式 为

$$f(x)=\sum _{i=0}^{n-1}\left(f\left(x_i\right)\prod _{\genfrac{}{}{0ex}{}{0\le j<n}{j\ne i}}\frac{x-x_j}{x_i-x_j}\right)$$

且 $x_i\ne x_j$ 对于所有 $i\ne j$ 成立。与 多项式的快速插值 中相同，我们令 $M(x):={\prod }_{i=0}^{n-1}\left(x-x_i\right)$，根据洛必达法则，有

$$M^{\prime }(x_i)=\underset{x\to x_i}{\lim }\frac{M(x)}{x-x_i}=\prod _{\genfrac{}{}{0ex}{}{0\le j<n}{j\ne i}}\left(x_i-x_j\right)$$

修正 Lagrange 插值公式 就是

$$f(x)=M(x)\left(\sum _{i=0}^{n-1}\frac{f\left(x_i\right)/M^{\prime }(x_i)}{x-x_i}\right)$$

那么现在我们有

$$f(x)=M(x)\left(\sum _{i=0}^{n-1}\frac{f\left(q^i\right)/M^{\prime }\left(q^i\right)}{x-q^i}\right)$$

其中 $M(x)={\prod }_{j=0}^{n-1}\left(x-q^j\right)$。若我们设 $n$ 为偶数，令 $n=2k$ 和 $H(x):={\prod }_{j=0}^{k-1}\left(x-q^j\right)$，那么

$$M(x)=H(x)\cdot q^{k^2}\cdot H\left(\frac{x}{q^k}\right)$$

这使得我们可以快速计算 $M(x)$。然后用 Bluestein 算法来计算 $M^{\prime }(1),\dots ,M^{\prime }(q^{n-1})$。令 $c_i:=f\left(q^i\right)/M^{\prime }\left(q^i\right)$，我们有

$$f(x)=M(x)\left(\sum _{i=0}^{n-1}\frac{c_i}{x-q^i}\right)$$

因为 $\mathrm{deg}f(x)<n$，我们只需计算 ${\sum }_{i=0}^{n-1}\frac{c_i}{x-q^i}\mathrm{mod}{x}^n$，其中 $\frac{c_i}{x-q^i}\in \mathbb{C}\left[\left[x\right]\right]$，也就是

$$\begin{array}{rl}\sum _{i=0}^{n-1}\frac{c_i}{x-q^i}\mathrm{mod}{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{array}$$

其中 $C(x)={\sum }_{i=0}^{n-1}{c}_i{x}^i$。我们可以用 Bluestein 算法来计算 $C\left(q^{-1}\right),\dots ,C\left(q^{-n}\right)$。

简单来说，我们分别进行下面的计算：

1. 用减治法（decrease and conquer）计算 $M(x)$；
2. 用 Bluestein 算法计算 $M^{\prime }(1),\dots ,M^{\prime }(q^{n-1})$；
3. 用 Bluestein 算法计算 $C\left(q^{-1}\right),\dots ,C\left(q^{-n}\right)$；
4. 用一次多项式乘法计算 $f(x)$。

其中每一步的时间复杂度都等于两个次数小于等于 $n$ 的多项式相乘的时间复杂度。

模板实现

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

参考文献

1. Bostan, A. (2010). Fast algorithms for polynomials and matrices. JNCF 2010. Algorithms Project, INRIA.