Biến đổi Chirp Z (CZT)

Tương tự biến đổi Fourier rời rạc, biến đổi Chirp Z là thuật toán cho đa thức \(f(x) = \sum_{i = 0}^{m - 1} f_i x^i \in \mathbb{C}\lbrack x\rbrack\) và \(q \in \mathbb{C} \setminus \{0\}\), tính \(f(1), f(q), \dots, f(q^{n - 1})\) mà không yêu cầu \(q\) là căn đơn vị. Cũng có thể dùng cho biến đổi số học. Phần sau sẽ giới thiệu Chirp Z và nghịch biến đổi.

Biến đổi Chirp Z

Theo định nghĩa, Chirp Z có thể viết:

\[ \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} \]

trong đó \(f(x) := \sum_{i = 0}^{m - 1} f_i x^i \in \mathbb{C}\lbrack x\rbrack\) và \(q \in \mathbb{C} \setminus \{0\}\).

Thuật toán Bluestein

Xét

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

với \(i, j \in \mathbb{Z}\), ta có thể xây dựng

\[ \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} \]

Trong đó \(G(x) \in \mathbb{C}\left\lbrack x, x^{-1}\right\rbrack\), và với \(i = 0, \dots, n - 1\) ta có

\[ \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} \]

và \(q^{\binom{i + 1}{2}} = q^{\binom{i}{2}}\cdot q^i\), \(\binom{-i}{2} = \binom{i + 1}{2}\). Như vậy chỉ cần một phép nhân đa thức, thuật toán này gọi là Bluestein.

Mẫu ( P6800【Mẫu】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;
}

Nghịch biến đổi Chirp Z

Nghịch Chirp Z có thể viết:

\[ \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) \]

trong đó \(f(x) \in \mathbb{C}\left\lbrack x\right\rbrack_{< n}\) và \(q \in \mathbb{C} \setminus \{0\}\), đồng thời \(q^i \neq q^j\) với mọi \(i \neq j\); đây là điều kiện nội suy đa thức.

Thuật toán Bostan–Schost

Nhắc lại công thức nội suy Lagrange:

\[ 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) \]

và \(x_i \neq x_j\) với mọi \(i \neq j\). Giống như trong nội suy đa thức nhanh, đặt \(M(x) := \prod_{i = 0}^{n - 1}\left(x - x_i\right)\), theo quy tắc L’Hôpital ta có

\[ 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) \]

Công thức Lagrange sửa đổi là

\[ f(x) = M(x)\left(\sum_{i = 0}^{n - 1}\frac{f\left(x_i\right)/M'(x_i)}{x - x_i}\right) \]

Do đó

\[ 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) \]

trong đó \(M(x)=\prod_{j = 0}^{n - 1}\left(x - q^j\right)\). Nếu giả sử \(n\) chẵn, đặt \(n = 2k\) và \(H(x) := \prod_{j = 0}^{k - 1}\left(x - q^j\right)\), thì

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

Nhờ đó ta có thể tính nhanh \(M(x)\). Sau đó dùng Bluestein để tính \(M'(1), \dots, M'(q^{n - 1})\). Đặt \(c_i := f\left(q^i\right)/M'\left(q^i\right)\), ta có

\[ f(x) = M(x)\left(\sum_{i = 0}^{n - 1}\frac{c_i}{x - q^i}\right) \]

Vì \(\deg f(x) < n\), ta chỉ cần tính \(\sum_{i = 0}^{n - 1}\frac{c_i}{x - q^i}\bmod{x^n}\), với \(\frac{c_i}{x - q^i} \in \mathbb{C}\left\lbrack\left\lbrack x\right\rbrack\right\rbrack\), tức là

\[ \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} \]

trong đó \(C(x) = \sum_{i = 0}^{n - 1} c_i x^i\). Ta có thể dùng Bluestein để tính \(C\left(q^{-1}\right), \dots, C\left(q^{-n}\right)\).

Tóm lại, ta thực hiện:

Dùng chia để trị (decrease and conquer) tính \(M(x)\);
Dùng Bluestein để tính \(M'(1), \dots, M'(q^{n - 1})\);
Dùng Bluestein để tính \(C\left(q^{-1}\right), \dots, C\left(q^{-n}\right)\);
Dùng một phép nhân đa thức để tính \(f(x)\).

Mỗi bước đều có độ phức tạp bằng thời gian nhân hai đa thức bậc không vượt quá \(n\).

Mẫu triển khai

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

Tài liệu tham khảo

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

Last updated on this page:, Update history
Found an error? Want to help improve? Edit this page on GitHub!
Contributors to this page:OI-wiki
All content on this page is provided under the terms of the CC BY-SA 4.0 and SATA license, additional terms may apply