Chirp Z Transform
(czt.hpp)

View this file on GitHub
Last update: 2025-01-19 15:28:01+08:00
Include: #include "czt.hpp"
Link: https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric

Chirp Z transform

The Chirp Z transform is by definition

\[\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}\left\lbrack x\right\rbrack$ and $q\in\mathbb{C}\setminus \left\lbrace 0\right\rbrace$. Since

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

note that $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 _ jq^{\binom{-j}{2} - \binom{i - j}{2}} \\ &= q^{-\binom{i}{2}}\cdot f\left(q^i\right) \end{aligned}\]

note that $q^{\binom{i+1}{2}}=q^{\binom{i}{2}}\cdot q^i$ and $\binom{-i}{2} = \binom{i+1}{2}$.

One can take advantage of the “middle product” algorithm. CZT could be done in time $O\left(\mathsf{M}(n+m)\right)$.

Inverse Chirp Z transform

The inverse Chirp Z transform is by definition

\[\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 _ {\lt n}$ and $q\in\mathbb{C}\setminus\left\lbrace 0\right\rbrace$, $q^{i}\neq q^{j}$ for all $i\neq j$.

(Modified) Lagrange interpolation formula

Recall the Lagrange interpolation formula, we have

\[f(x) = \sum _ {i=0}^{n-1}\left(f\left(x _ i\right)\prod _ {0\leq j\lt n\atop j\neq i} \frac{x - x _ j}{x _ i - x _ j}\right)\]

for $x _ i \neq x _ j$ for all $i\neq j$.

L’Hôpital’s Rule: The derivative of $f = \sum _ {i} f _ i x^i$ is $f’ = \sum _ {i} i f_i x^{i - 1}$, for $g\in\mathbb{C}\left\lbrack x\right\rbrack$ we have the product rule: $\left(fg\right)’ = f’g + fg’$. Let $\alpha\in\mathbb{C}$ and $f(\alpha) = 0$ then $\left(fg\right)’(\alpha) = f’(\alpha)g(\alpha)$, if $f’(\alpha)\neq 0$ then $(fg/f)(\alpha) = g(\alpha) = (fg)’(\alpha)/f’(\alpha)$.

Let $M(x):=\prod _ {i=0}^{n-1}\left(x - x _ i\right)$, we have

\[M'(x _ i)=\prod_{0\leq j\lt n\atop j\neq i}\left(x _ i - x _ j\right)\]

The modified Lagrange interpolation formula is

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

Now we have

\[f(x)=M(x)\cdot\sum _ {i=0}^{n-1}\frac{f\left(q^i\right)/M'\left(q^i\right)}{x-q^i}\]

where $M(x)=\prod _ {j=0}^{n-1}\left(x - q^j\right)$. If we set $n-1=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)\]

So we could compute $M(x)$ in time $O\left(\mathsf{M}(n)\right)$. We could use CZT 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)\cdot\sum _ {i=0}^{n-1}\frac{c _ i}{x-q^i}\]

since $\deg f(x) \lt n$, we only need to compute $\sum _ {i=0}^{n-1}\frac{c _ i}{x - q^i}\bmod{x^n}$, which 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 another CZT to compute $C\left(q^{-1}\right),\dots ,C\left(q^{-n}\right)$.

// An alternative implementation of this file.
template<typename Tp> inline std::vector<Tp> inv_czt(const std::vector<Tp> &f, const Tp q) {
    if (f.empty()) return {};
    const int n = f.size();
    // prod_(i=0..(n-1))(x-q^i)
    auto rec = [q, iq = q.inv()](auto &&rec, int n) -> std::vector<Tp> {
        if (n == 0) return std::vector<Tp>{Tp(1)};
        if (n == 1) return std::vector<Tp>{Tp(-1), Tp(1)};
        const auto H   = rec(rec, n / 2);
        auto HH        = H;
        const auto qn  = q.pow(n / 2);
        const auto iqn = qn.inv();
        Tp qq          = q.pow((long long)(n / 2) * (n / 2));
        for (int i = 0; i <= n / 2; ++i) HH[i] *= qq, qq *= iqn;
        const auto res = convolution(H, HH);
        return (n & 1) ? convolution(res, std::vector<Tp>{-q.pow(n - 1), Tp(1)}) : res;
    };

    const auto M = rec(rec, n);
    auto C       = czt(deriv(M), q, n);
    for (int i = 0; i < n; ++i) C[i] = f[i] / C[i];
    C = czt(C, q.inv(), n, q.inv());
    for (int i = 0; i < n; ++i) C[i] = -C[i];
    auto res = convolution(M, C);
    res.resize(n);
    return res;
}

noshi91 showed that we could compute $M(x)$ and $M’(1),\dots ,M’(q^{n-1})$ faster. Let $s _ i:=\prod _ {j=1}^i\left(1-q^j\right)$ and $s _ 0 := 0$ followed by noshi91’s definition. We have

\[\begin{aligned} \prod _ {0\leq j\lt n\atop j\neq i}\left(q^i-q^j\right)&=\left(\prod _ {j=0}^{i-1}\left(q^i-q^j\right)\right)\left(\prod _ {j=i+1}^{n-1}\left(q^i-q^j\right)\right) \\ &=\left(\prod _ {j=0}^{i-1}q^j\left(q^{i-j}-1\right)\right)\left(\prod _ {j=i+1}^{n-1}q^i\left(1-q^{j-i}\right)\right) \\ &=(-1)^iq^{\sum _ {j=0}^{i-1}j}\left(\prod _ {j=0}^{i-1}\left(1-q^{i-j}\right)\right)\cdot q^{i(n-i-1)}\left(\prod _ {j=i+1}^{n-1}\left(1-q^{j-i}\right)\right) \\ &=(-1)^iq^{\binom{i}{2}}\left(\prod _ {k=1}^{i}\left(1-q^k\right)\right)\cdot q^{i(n-i-1)}\left(\prod _ {k=1}^{n-i-1}\left(1-q^k\right)\right) \\ &=(-1)^iq^{\binom{i}{2}}s _ i\cdot q^{i(n-i-1)}s _ {n-i-1} \end{aligned}\]

note that $q^{\binom{i+1}{2}}\cdot q^{(i+1)(n-(i+1)-1)}=\left(q^{\binom{i}{2}}\cdot q^{i}\right)\left(q^{i(n-i-1)}\cdot q^{n-2i-2}\right)$. Now we only need to compute $M(x)$.

$q$-analog

Let $n\in\mathbb{N}$, the $q$-analog of $n$ is defined by

\[\left\lbrack n\right\rbrack _ q := \begin{cases} 0,&\text{if }n=0, \\ 1+q+\cdots +q^{n-1},&\text{otherwise}. \end{cases}\]

note that $\left\lbrack n\right\rbrack _ q = n$ when $q = 1$. And $q\left\lbrack n\right\rbrack _ q =q+q^2+\cdots +q^n$, so

\[\left\lbrack n\right\rbrack _ q= \begin{cases} n,&\text{if }q=1, \\ \frac{1-q^n}{1-q},&\text{otherwise}. \end{cases}\]

We could also findout that

\[\begin{aligned} \left\lbrack -n\right\rbrack _ q &= \frac{1-q^{-n}}{1-q} \\ &= q^{-n}\frac{q^n-1}{1-q} \\ &= -q^{-n}\left\lbrack n\right\rbrack _ q \end{aligned}\]

Now we define the $q$-factorial

\[n! _ q := \begin{cases} 1,&\text{if }n=0, \\ \left\lbrack 1\right\rbrack_q\left\lbrack 2\right\rbrack _ q\cdots \left\lbrack n\right\rbrack _ q,&\text{otherwise}. \end{cases}\]

Finally we could define the $q$-binomial coefficients

\[\binom{n}{k} _ q := \begin{cases} \frac{n! _ q}{k! _ q(n-k)! _ q},&\text{if }0\leq k\leq n, \\ 0,&\text{otherwise}. \end{cases}\]

Before introducing the $q$-binomial theorem, let’s show some common identities. If $q \neq 1$, we have

\[\begin{aligned} \left\lbrack n\right\rbrack _ q&=\frac{1-q^n}{1-q} \\ &=\frac{1-q^k+q^k-q^n}{1-q} \\ &=\frac{1-q^k}{1-q}+q^k\frac{1-q^{n-k}}{1-q} \\ &=\left\lbrack k\right\rbrack _ q +q^k\left\lbrack n-k\right\rbrack _ q \end{aligned}\]

which still works when $q = 1$ since $n = k + (n - k)$. We have

\[\begin{aligned} \binom{n+1}{k} _ q&=\frac{(n+1)! _ q}{k! _ q(n+1-k)! _ q} \\ &=\frac{n! _ q}{k! _ q(n+1-k)! _ q}\cdot \left(\left\lbrack k\right\rbrack _ q +q^k\left\lbrack n+1-k\right\rbrack _ q\right) \\ &=\frac{n! _ q}{(k-1)! _ q(n+1-k)! _ q}\cdot \left\lbrack k\right\rbrack _ q^{-1}\cdot \left(\left\lbrack k\right\rbrack _ q +q^k\left\lbrack n+1-k\right\rbrack _ q\right) \\ &=\frac{n! _ q}{(k-1)! _ q(n+1-k)! _ q}+q^k\frac{n! _ q\left\lbrack n+1-k\right\rbrack _ q}{k! _ q(n+1-k)! _ q} \\ &=\binom{n}{k-1} _ q+q^k\binom{n}{k} _ q \end{aligned}\]

Since $\binom{n}{k} _ q=\binom{n}{n-k} _ q$, we could replace $k$ by $n-k$ which gives

\[\begin{aligned} \binom{n+1}{k} _ q&=\binom{n+1}{n+1-k} _ q \\ &=\binom{n}{n-k} _ q+q^{n+1-k}\binom{n}{n+1-k} _ q \\ &=\binom{n}{k} _ q+q^{n-k+1}\binom{n}{k-1} _ q \end{aligned}\]

These are called the $q$-Pascal relation.

Rothe’s $q$-binomial Theorem

For variable $q,a,x$, we have

\[\prod _ {i=0}^{n-1}\left(a+q^ix\right) = \sum _ {k=0}^n\binom{n}{k} _ q q^{\binom{k}{2}}a^{n-k}x^k\]

and LHS is defined to be $1$ when $n = 0$.

Proof. Let RHS equals $r_n(x,a)$, apply $q$-Pascal relation here, we have

\[\begin{aligned} r _ {n+1}(x,a)&=\sum _ {k=0}^{n+1}\binom{n+1}{k} _ q q^{\binom{k}{2}}x^ka^{n+1-k} \\ &=\sum _ {k=0}^{n+1}\left(\binom{n}{k-1} _ q+q^k\binom{n}{k} _ q\right)q^{\binom{k}{2}}x^ka^{n+1-k} \\ &=\left(\sum _ {k=1}^{n+1}\binom{n}{k-1} _ q q^{\binom{k}{2}}x^ka^{n+1-k}\right)+\left(\sum _ {k=0}^{n}\binom{n}{k} _ q q^{\binom{k}{2}}\left(qx\right)^ka^{n+1-k}\right) \\ &=\left(\sum _ {j=0}^{n}\binom{n}{j} _ qq^{\binom{j+1}{2}}x^{j+1}a^{n-j}\right)+\left(\sum _ {j=0}^{n}\binom{n}{j} _ q q^{\binom{j}{2}}\left(qx\right)^ja^{n+1-j}\right) \end{aligned}\]

Since $\binom{n}{-1} _ q = \binom{n}{n+1} _ q = 0$ and $\binom{j+1}{2}=\binom{j}{2}+j$ for $j\in\mathbb{N}$, we have

\[\begin{aligned} r _ {n+1}(x,a)&=\left(\sum _ {j=0}^{n}\binom{n}{j} _ q q^{\binom{j}{2}}q^jx^{j+1}a^{n-j}\right)+\left(\sum _ {j=0}^{n}\binom{n}{j} _ q q^{\binom{j}{2}}\left(qx\right)^ja^{n+1-j}\right) \\ &=\left(\sum _ {j=0}^{n}\binom{n}{j} _ q q^{\binom{j}{2}}\left(qx\right)^ja^{n-j}x\right)+\left(\sum _ {j=0}^{n}\binom{n}{j} _ q q^{\binom{j}{2}}\left(qx\right)^ja^{n+1-j}\right) \\ &=\sum _ {j=0}^{n}\binom{n}{j} _ q q^{\binom{j}{2}}\left(qx\right)^ja^{n-j}(x+a) \\ &=(a+x)r _ n(qx,a) \end{aligned}\]

By induction, we have

\[\begin{aligned} r _ {n+1}(x,a)&=(a+x)r _ n(qx,a) \\ &=(a+x)(a+qx)r_{n-1}\left(q^2x,a\right) \\ &=\cdots \\ &=(a+x)(a+qx)\cdots \left(a+q^nx\right)r _ 0\left(q^{n+1}x,a\right) \end{aligned}\]

and $r_0(u,v) = 1$ agrees with our definition.

Now back to our problem: Computing $M(x)$ mentioned above.

\[\begin{aligned} M(x)&=\prod _ {j=0}^{n-1}\left(x-q^j\right) \\ &=\sum _ {k=0}^n\binom{n}{k} _ q q^{\binom{k}{2}}(-1)^kx^{n-k} \end{aligned}\]

where

\[\binom{n}{k} _ q=\frac{\left(1-q\right)\cdots \left(1-q^n\right)}{\left(1-q\right)\cdots \left(1-q^k\right)\cdot \left(1-q\right)\cdots \left(1-q^{n-k}\right)}=\frac{s _ n}{s _ k s _ {n-k}}\]

Since $q^n \neq 1$ is not guaranteed, $M(0)$ should be computed separately.

References

37zigen. 多項式補間：アルゴリズム. url: https://37zigen.com/lagrange-interpolation/
noshi91. 標本点が等比数列を成す場合に補間多項式を計算するアルゴリズム. url: https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric
Bostan, A. (2010). Fast algorithms for polynomials and matrices. JNCF 2010. Algorithms Project, INRIA. url: https://specfun.inria.fr/bostan/publications/exposeJNCF.pdf
Warren P. Johnson. An Introduction to $q$-analysis. American Mathematical Soc., 2020.

Depends on

Verified with

Code

#pragma once

#include "batch_inv.hpp"
#include "middle_product.hpp"
#include "poly_basic.hpp"
#include <algorithm>
#include <vector>

// returns F(a),F(ac),F(ac^2),...,F(ac^(n-1))
// Use        ij = binom(i,2)   + binom(-j,2) - binom(i-j,2)
// instead of ij = binom(i+j,2) - binom(i,2)  - binom(j,2)
template<typename Tp> inline std::vector<Tp> czt(std::vector<Tp> F, Tp c, int n, Tp a = 1) {
    if (n <= 0) return {};
    const int degF = degree(F);
    shrink(F);
    if (degF < 0) return std::vector<Tp>(n);
    if (degF == 0 || a == 0) return std::vector<Tp>(n, F[0]);
    if (a != 1) {
        // F(x) <- F(ax)
        Tp aa = 1;
        for (int i = 0; i <= degF; ++i) F[i] *= aa, aa *= a;
    }
    if (c == 0) {
        std::vector<Tp> res(n, F[0]);
        for (int i = 1; i <= degF; ++i) res[0] += F[i];
        return res;
    }
    std::vector<Tp> H(std::max(degF + 1, n - 1));
    Tp cc = H[0] = 1;
    for (int i = 1; i < (int)H.size(); ++i) H[i] = H[i - 1] * (cc *= c);
    std::vector<Tp> G(degF + n); // G[i+degF]=c^(-binom(i,2))
    auto GG     = G.begin() + degF;
    const Tp ic = c.inv();
    cc = GG[0] = 1;
    for (int i = 1; i < n; ++i) GG[i] = GG[i - 1] * cc, cc *= ic;
    cc = 1;
    for (int i = -1; i >= -degF; --i) GG[i] = GG[i + 1] * (cc *= ic);
    // F[i] <- c^(binom(i+1,2))*F[i]
    for (int i = 0; i <= degF; ++i) F[i] *= H[i];
    F = middle_product(G, F);
    // F[i] <- c^(binom(i,2))*F[i]
    for (int i = 1; i < n; ++i) F[i] *= H[i - 1];
    return F;
}

// returns f s.t. f(aq^i)=F[i]
// aq^i != aq^j for all i != j
// see: https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric
// noshi91. 標本点が等比数列を成す場合に補間多項式を計算するアルゴリズム.
template<typename Tp> inline std::vector<Tp> inv_czt(const std::vector<Tp> &F, Tp q, Tp a = 1) {
    if (F.empty()) return {};
    if (a == 0) return {F[0]};
    const int n = F.size();
    std::vector<Tp> Q(n), S(n), M(n), D(n);
    Tp qq = 1;
    // Q[i]=q^i
    for (int i = 0; i < n; ++i) Q[i] = qq, qq *= q;
    // S[i]=prod_(i=1..i)(1-q^i)
    S[0] = 1;
    for (int i = 1; i < n; ++i) S[i] = S[i - 1] * (1 - Q[i]);
    const auto Sn   = S[n - 1] * (1 - qq);
    const auto invS = batch_inv(S);
    qq              = 1;
    // M[n-i]=qbinom(n,i)*q^(binom(i,2))*(-1)^i
    for (int i = 1; i < n; ++i) M[n - i] = Sn * invS[i] * invS[n - i] * (qq *= -Q[i - 1]);
    M[0] = qq * -Q[n - 1]; // in case of q^n=1
    // D[i]=S[i]*S[n-i-1]*q^(binom(i,2)+i(n-i-1))*(-1)^i
    D[0] = 1;
    for (int i = 0; i < n - 1; ++i) D[i + 1] = D[i] * -Q[n - i - 2];
    for (int i = 0; i < n; ++i) D[i] *= S[i] * S[n - i - 1];
    // D[i] <- -F[i]/D[i]
    D = batch_inv(D);
    for (int i = 0; i < n; ++i) D[i] *= -F[i];
    auto res = convolution(M, czt(D, q.inv(), n, q.inv()));
    res.resize(n);
    if (a != 1) {
        const auto ia = a.inv();
        Tp aa         = 1;
        for (int i = 0; i < n; ++i) res[i] *= aa, aa *= ia;
    }
    return res;
}

#line 2 "czt.hpp"

#line 2 "batch_inv.hpp"

#include <cassert>
#include <vector>

template<typename Tp> inline std::vector<Tp> batch_inv(const std::vector<Tp> &a) {
    if (a.empty()) return {};
    const int n = a.size();
    std::vector<Tp> b(n);
    Tp v = 1;
    for (int i = 0; i < n; ++i) b[i] = v, v *= a[i];
    assert(v != 0);
    v = v.inv();
    for (int i = n - 1; i >= 0; --i) b[i] *= v, v *= a[i];
    return b;
}
#line 2 "middle_product.hpp"

#line 2 "fft.hpp"

#include <algorithm>
#line 5 "fft.hpp"
#include <iterator>
#include <memory>
#line 8 "fft.hpp"

template<typename Tp> class FftInfo {
    static Tp least_quadratic_nonresidue() {
        for (int i = 2;; ++i)
            if (Tp(i).pow((Tp::mod() - 1) / 2) == -1) return Tp(i);
    }

    const int ordlog2_;
    const Tp zeta_;
    const Tp invzeta_;
    const Tp imag_;
    const Tp invimag_;

    mutable std::vector<Tp> root_;
    mutable std::vector<Tp> invroot_;

    FftInfo()
        : ordlog2_(__builtin_ctzll(Tp::mod() - 1)),
          zeta_(least_quadratic_nonresidue().pow((Tp::mod() - 1) >> ordlog2_)),
          invzeta_(zeta_.inv()), imag_(zeta_.pow(1LL << (ordlog2_ - 2))), invimag_(-imag_),
          root_{Tp(1), imag_}, invroot_{Tp(1), invimag_} {}

public:
    static const FftInfo &get() {
        static FftInfo info;
        return info;
    }

    Tp imag() const { return imag_; }
    Tp inv_imag() const { return invimag_; }
    Tp zeta() const { return zeta_; }
    Tp inv_zeta() const { return invzeta_; }
    const std::vector<Tp> &root(int n) const {
        // [0, n)
        assert((n & (n - 1)) == 0);
        if (const int s = root_.size(); s < n) {
            root_.resize(n);
            for (int i = __builtin_ctz(s); (1 << i) < n; ++i) {
                const int j = 1 << i;
                root_[j]    = zeta_.pow(1LL << (ordlog2_ - i - 2));
                for (int k = j + 1; k < j * 2; ++k) root_[k] = root_[k - j] * root_[j];
            }
        }
        return root_;
    }
    const std::vector<Tp> &inv_root(int n) const {
        // [0, n)
        assert((n & (n - 1)) == 0);
        if (const int s = invroot_.size(); s < n) {
            invroot_.resize(n);
            for (int i = __builtin_ctz(s); (1 << i) < n; ++i) {
                const int j = 1 << i;
                invroot_[j] = invzeta_.pow(1LL << (ordlog2_ - i - 2));
                for (int k = j + 1; k < j * 2; ++k) invroot_[k] = invroot_[k - j] * invroot_[j];
            }
        }
        return invroot_;
    }
};

inline int fft_len(int n) {
    --n;
    n |= n >> 1, n |= n >> 2, n |= n >> 4, n |= n >> 8;
    return (n | n >> 16) + 1;
}

namespace detail {

template<typename Iterator> inline void
butterfly_n(Iterator a, int n,
            const std::vector<typename std::iterator_traits<Iterator>::value_type> &root) {
    assert(n > 0);
    assert((n & (n - 1)) == 0);
    const int bn = __builtin_ctz(n);
    if (bn & 1) {
        for (int i = 0; i < n / 2; ++i) {
            const auto a0 = a[i], a1 = a[i + n / 2];
            a[i] = a0 + a1, a[i + n / 2] = a0 - a1;
        }
    }
    for (int i = n >> (bn & 1); i >= 4; i /= 4) {
        const int i4 = i / 4;
        for (int k = 0; k < i4; ++k) {
            const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
            const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
            const auto a02p = a0 + a2, a02m = a0 - a2;
            const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1];
            a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p;
            a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m;
        }
        for (int j = i, m = 2; j < n; j += i, m += 2) {
            const auto r = root[m], r2 = r * r, r3 = r2 * r;
            for (int k = j; k < j + i4; ++k) {
                const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1] * r;
                const auto a2 = a[k + i4 * 2] * r2, a3 = a[k + i4 * 3] * r3;
                const auto a02p = a0 + a2, a02m = a0 - a2;
                const auto a13p = a1 + a3, a13m = (a1 - a3) * root[1];
                a[k + i4 * 0] = a02p + a13p, a[k + i4 * 1] = a02p - a13p;
                a[k + i4 * 2] = a02m + a13m, a[k + i4 * 3] = a02m - a13m;
            }
        }
    }
}

template<typename Iterator> inline void
inv_butterfly_n(Iterator a, int n,
                const std::vector<typename std::iterator_traits<Iterator>::value_type> &root) {
    assert(n > 0);
    assert((n & (n - 1)) == 0);
    const int bn = __builtin_ctz(n);
    for (int i = 4; i <= (n >> (bn & 1)); i *= 4) {
        const int i4 = i / 4;
        for (int k = 0; k < i4; ++k) {
            const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
            const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
            const auto a01p = a0 + a1, a01m = a0 - a1;
            const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1];
            a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = a01m + a23m;
            a[k + i4 * 2] = a01p - a23p, a[k + i4 * 3] = a01m - a23m;
        }
        for (int j = i, m = 2; j < n; j += i, m += 2) {
            const auto r = root[m], r2 = r * r, r3 = r2 * r;
            for (int k = j; k < j + i4; ++k) {
                const auto a0 = a[k + i4 * 0], a1 = a[k + i4 * 1];
                const auto a2 = a[k + i4 * 2], a3 = a[k + i4 * 3];
                const auto a01p = a0 + a1, a01m = a0 - a1;
                const auto a23p = a2 + a3, a23m = (a2 - a3) * root[1];
                a[k + i4 * 0] = a01p + a23p, a[k + i4 * 1] = (a01m + a23m) * r;
                a[k + i4 * 2] = (a01p - a23p) * r2, a[k + i4 * 3] = (a01m - a23m) * r3;
            }
        }
    }
    if (bn & 1) {
        for (int i = 0; i < n / 2; ++i) {
            const auto a0 = a[i], a1 = a[i + n / 2];
            a[i] = a0 + a1, a[i + n / 2] = a0 - a1;
        }
    }
}

} // namespace detail

// FFT_n: A(x) |-> bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))]
template<typename Iterator> inline void fft_n(Iterator a, int n) {
    using Tp = typename std::iterator_traits<Iterator>::value_type;
    detail::butterfly_n(a, n, FftInfo<Tp>::get().root(n / 2));
}

template<typename Tp> inline void fft(std::vector<Tp> &a) { fft_n(a.begin(), a.size()); }

// IFFT_n: bit-reversed order of [A(1), A(zeta_n), ..., A(zeta_n^(n-1))] |-> A(x)
template<typename Iterator> inline void inv_fft_n(Iterator a, int n) {
    using Tp = typename std::iterator_traits<Iterator>::value_type;
    detail::inv_butterfly_n(a, n, FftInfo<Tp>::get().inv_root(n / 2));
    const Tp iv = Tp::mod() - (Tp::mod() - 1) / n;
    for (int i = 0; i < n; ++i) a[i] *= iv;
}

template<typename Tp> inline void inv_fft(std::vector<Tp> &a) { inv_fft_n(a.begin(), a.size()); }

// IFFT_n^T: A(x) |-> 1/n FFT_n((x^n A(x^(-1))) mod (x^n - 1))
template<typename Iterator> inline void transposed_inv_fft_n(Iterator a, int n) {
    using Tp    = typename std::iterator_traits<Iterator>::value_type;
    const Tp iv = Tp::mod() - (Tp::mod() - 1) / n;
    for (int i = 0; i < n; ++i) a[i] *= iv;
    detail::butterfly_n(a, n, FftInfo<Tp>::get().inv_root(n / 2));
}

template<typename Tp> inline void transposed_inv_fft(std::vector<Tp> &a) {
    transposed_inv_fft_n(a.begin(), a.size());
}

// FFT_n^T : FFT_n((x^n A(x^(-1))) mod (x^n - 1)) |-> n A(x)
template<typename Iterator> inline void transposed_fft_n(Iterator a, int n) {
    using Tp = typename std::iterator_traits<Iterator>::value_type;
    detail::inv_butterfly_n(a, n, FftInfo<Tp>::get().root(n / 2));
}

template<typename Tp> inline void transposed_fft(std::vector<Tp> &a) {
    transposed_fft_n(a.begin(), a.size());
}

template<typename Tp> inline std::vector<Tp> convolution_fft(std::vector<Tp> a, std::vector<Tp> b) {
    if (a.empty() || b.empty()) return {};
    const int n   = a.size();
    const int m   = b.size();
    const int len = fft_len(n + m - 1);
    a.resize(len);
    b.resize(len);
    fft(a);
    fft(b);
    for (int i = 0; i < len; ++i) a[i] *= b[i];
    inv_fft(a);
    a.resize(n + m - 1);
    return a;
}

template<typename Tp> inline std::vector<Tp> square_fft(std::vector<Tp> a) {
    if (a.empty()) return {};
    const int n   = a.size();
    const int len = fft_len(n * 2 - 1);
    a.resize(len);
    fft(a);
    for (int i = 0; i < len; ++i) a[i] *= a[i];
    inv_fft(a);
    a.resize(n * 2 - 1);
    return a;
}

template<typename Tp>
inline std::vector<Tp> convolution_naive(const std::vector<Tp> &a, const std::vector<Tp> &b) {
    if (a.empty() || b.empty()) return {};
    const int n = a.size();
    const int m = b.size();
    std::vector<Tp> res(n + m - 1);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j) res[i + j] += a[i] * b[j];
    return res;
}

template<typename Tp>
inline std::vector<Tp> convolution(const std::vector<Tp> &a, const std::vector<Tp> &b) {
    if (std::min(a.size(), b.size()) < 60) return convolution_naive(a, b);
    if (std::addressof(a) == std::addressof(b)) return square_fft(a);
    return convolution_fft(a, b);
}
#line 7 "middle_product.hpp"

// see:
// [1]: Guillaume Hanrot, Michel Quercia, Paul Zimmermann. The Middle Product Algorithm I.
// [2]: Alin Bostan, Grégoire Lecerf, Éric Schost. Tellegen's principle into practice.

template<typename Tp>
inline std::vector<Tp> middle_product_naive(const std::vector<Tp> &f, const std::vector<Tp> &g) {
    const int m = f.size();
    const int n = g.size();
    assert(m >= n);
    std::vector<Tp> res(m - n + 1);
    for (int i = n - 1; i < m; ++i)
        for (int j = i - (n - 1); j <= i; ++j) res[i - (n - 1)] += f[j] * g[i - j];
    return res;
}

template<typename Tp>
inline std::vector<Tp> middle_product_fft(std::vector<Tp> f, std::vector<Tp> g) {
    const int m = f.size();
    const int n = g.size();
    assert(m >= n);
    std::reverse(g.begin(), g.end());
    const int len = fft_len(m);
    f.resize(len);
    g.resize(len);
    transposed_inv_fft(f);
    fft(g);
    for (int i = 0; i < len; ++i) f[i] *= g[i];
    transposed_fft(f);
    f.resize(m - n + 1);
    return f;
}

// returns (fg)_(n-1),...,(fg)_(m-1)
// f: f_0 + ... + f_(m-1)x^(m-1)
// g: g_0 + ... + g_(n-1)x^(n-1)
// requires m >= n
template<typename Tp>
inline std::vector<Tp> middle_product(const std::vector<Tp> &f, const std::vector<Tp> &g) {
    assert(f.size() >= g.size());
    if (f.size() < 60) return middle_product_naive(f, g);
    return middle_product_fft(f, g);
}
#line 2 "poly_basic.hpp"

#line 2 "binomial.hpp"

#line 5 "binomial.hpp"

template<typename Tp> class Binomial {
    std::vector<Tp> factorial_, invfactorial_;

    Binomial() : factorial_{Tp(1)}, invfactorial_{Tp(1)} {}

    void preprocess(int n) {
        if (const int nn = factorial_.size(); nn < n) {
            int k = nn;
            while (k < n) k *= 2;
            k = std::min<long long>(k, Tp::mod());
            factorial_.resize(k);
            invfactorial_.resize(k);
            for (int i = nn; i < k; ++i) factorial_[i] = factorial_[i - 1] * i;
            invfactorial_.back() = factorial_.back().inv();
            for (int i = k - 2; i >= nn; --i) invfactorial_[i] = invfactorial_[i + 1] * (i + 1);
        }
    }

public:
    static const Binomial &get(int n) {
        static Binomial bin;
        bin.preprocess(n);
        return bin;
    }

    Tp binom(int n, int m) const {
        return n < m ? Tp() : factorial_[n] * invfactorial_[m] * invfactorial_[n - m];
    }
    Tp inv(int n) const { return factorial_[n - 1] * invfactorial_[n]; }
    Tp factorial(int n) const { return factorial_[n]; }
    Tp inv_factorial(int n) const { return invfactorial_[n]; }
};
#line 2 "fps_basic.hpp"

#line 2 "semi_relaxed_conv.hpp"

#line 5 "semi_relaxed_conv.hpp"
#include <type_traits>
#include <utility>
#line 8 "semi_relaxed_conv.hpp"

template<typename Tp, typename Closure>
inline std::enable_if_t<std::is_invocable_r_v<Tp, Closure, int, const std::vector<Tp> &>,
                        std::vector<Tp>>
semi_relaxed_convolution_naive(const std::vector<Tp> &A, Closure gen, int n) {
    std::vector<Tp> B(n), AB(n);
    for (int i = 0; i < n; ++i) {
        for (int j = std::max(0, i - (int)A.size() + 1); j < i; ++j) AB[i] += A[i - j] * B[j];
        B[i] = gen(i, AB);
        if (!A.empty()) AB[i] += A[0] * B[i];
    }
    return B;
}

// returns coefficients generated by closure
// closure: gen(index, current_product)
template<typename Tp, typename Closure>
inline std::enable_if_t<std::is_invocable_r_v<Tp, Closure, int, const std::vector<Tp> &>,
                        std::vector<Tp>>
semi_relaxed_convolution(const std::vector<Tp> &A, Closure gen, int n) {
    if (A.size() < 60) return semi_relaxed_convolution_naive(A, gen, n);
    enum { BaseCaseSize = 32 };
    static_assert((BaseCaseSize & (BaseCaseSize - 1)) == 0);

    static const int Block[]     = {16, 16, 16, 16, 16};
    static const int BlockSize[] = {
        BaseCaseSize,
        BaseCaseSize * Block[0],
        BaseCaseSize * Block[0] * Block[1],
        BaseCaseSize * Block[0] * Block[1] * Block[2],
        BaseCaseSize * Block[0] * Block[1] * Block[2] * Block[3],
        BaseCaseSize * Block[0] * Block[1] * Block[2] * Block[3] * Block[4],
    };

    // returns (which_block, level)
    auto blockinfo = [](int ind) {
        int i = ind / BaseCaseSize, lv = 0;
        while ((i & (Block[lv] - 1)) == 0) i /= Block[lv++];
        return std::make_pair(i & (Block[lv] - 1), lv);
    };

    std::vector<Tp> B(n), AB(n);
    std::vector<std::vector<std::vector<Tp>>> dftA, dftB;

    for (int i = 0; i < n; ++i) {
        const int s = i & (BaseCaseSize - 1);

        // block contribution
        if (i >= BaseCaseSize && s == 0) {
            const auto [j, lv]  = blockinfo(i);
            const int blocksize = BlockSize[lv];

            if (blocksize * j == i) {
                if ((int)dftA.size() == lv) {
                    dftA.emplace_back();
                    dftB.emplace_back(Block[lv] - 1);
                }
                if ((j - 1) * blocksize < (int)A.size()) {
                    dftA[lv]
                        .emplace_back(A.begin() + (j - 1) * blocksize,
                                      A.begin() + std::min<int>((j + 1) * blocksize, A.size()))
                        .resize(blocksize * 2);
                    fft(dftA[lv][j - 1]);
                }
            }

            if (!dftA[lv].empty()) {
                dftB[lv][j - 1].resize(blocksize * 2);
                std::copy_n(B.begin() + (i - blocksize), blocksize, dftB[lv][j - 1].begin());
                std::fill_n(dftB[lv][j - 1].begin() + blocksize, blocksize, Tp(0));
                fft(dftB[lv][j - 1]);

                // middle product
                std::vector<Tp> mp(blocksize * 2);
                for (int k = 0; k < std::min<int>(j, dftA[lv].size()); ++k)
                    for (int l = 0; l < blocksize * 2; ++l)
                        mp[l] += dftA[lv][k][l] * dftB[lv][j - 1 - k][l];
                inv_fft(mp);

                for (int k = 0; k < blocksize && i + k < n; ++k) AB[i + k] += mp[k + blocksize];
            }
        }

        // basecase contribution
        for (int j = std::max(i - s, i - (int)A.size() + 1); j < i; ++j) AB[i] += A[i - j] * B[j];
        B[i] = gen(i, AB);
        if (!A.empty()) AB[i] += A[0] * B[i];
    }

    return B;
}
#line 8 "fps_basic.hpp"

template<typename Tp> inline int order(const std::vector<Tp> &a) {
    for (int i = 0; i < (int)a.size(); ++i)
        if (a[i] != 0) return i;
    return -1;
}

template<typename Tp> inline std::vector<Tp> fps_inv(const std::vector<Tp> &a, int n) {
    assert(order(a) == 0);
    if (n <= 0) return {};
    if (std::min<int>(a.size(), n) < 60)
        return semi_relaxed_convolution(
            a, [v = a[0].inv()](int n, auto &&c) { return n == 0 ? v : -c[n] * v; }, n);
    enum { Threshold = 32 };
    const int len = fft_len(n);
    std::vector<Tp> invA, shopA(len), shopB(len);
    invA = semi_relaxed_convolution(
        a, [v = a[0].inv()](int n, auto &&c) { return n == 0 ? v : -c[n] * v; }, Threshold);
    invA.resize(len);
    for (int i = Threshold * 2; i <= len; i *= 2) {
        std::fill(std::copy_n(a.begin(), std::min<int>(a.size(), i), shopA.begin()),
                  shopA.begin() + i, Tp(0));
        std::copy_n(invA.begin(), i, shopB.begin());
        fft_n(shopA.begin(), i);
        fft_n(shopB.begin(), i);
        for (int j = 0; j < i; ++j) shopA[j] *= shopB[j];
        inv_fft_n(shopA.begin(), i);
        std::fill_n(shopA.begin(), i / 2, Tp(0));
        fft_n(shopA.begin(), i);
        for (int j = 0; j < i; ++j) shopA[j] *= shopB[j];
        inv_fft_n(shopA.begin(), i);
        for (int j = i / 2; j < i; ++j) invA[j] = -shopA[j];
    }
    invA.resize(n);
    return invA;
}

template<typename Tp>
inline std::vector<Tp> fps_div(const std::vector<Tp> &a, const std::vector<Tp> &b, int n) {
    assert(order(b) == 0);
    if (n <= 0) return {};
    return semi_relaxed_convolution(
        b,
        [&, v = b[0].inv()](int n, auto &&c) {
            if (n < (int)a.size()) return (a[n] - c[n]) * v;
            return -c[n] * v;
        },
        n);
}

template<typename Tp> inline std::vector<Tp> deriv(const std::vector<Tp> &a) {
    const int n = (int)a.size() - 1;
    if (n <= 0) return {};
    std::vector<Tp> res(n);
    for (int i = 1; i <= n; ++i) res[i - 1] = a[i] * i;
    return res;
}

template<typename Tp> inline std::vector<Tp> integr(const std::vector<Tp> &a, Tp c = {}) {
    const int n = a.size() + 1;
    auto &&bin  = Binomial<Tp>::get(n);
    std::vector<Tp> res(n);
    res[0] = c;
    for (int i = 1; i < n; ++i) res[i] = a[i - 1] * bin.inv(i);
    return res;
}

template<typename Tp> inline std::vector<Tp> fps_log(const std::vector<Tp> &a, int n) {
    return integr(fps_div(deriv(a), a, n - 1));
}

template<typename Tp> inline std::vector<Tp> fps_exp(const std::vector<Tp> &a, int n) {
    if (n <= 0) return {};
    assert(a.empty() || a[0] == 0);
    return semi_relaxed_convolution(
        deriv(a),
        [bin = Binomial<Tp>::get(n)](int n, auto &&c) {
            return n == 0 ? Tp(1) : c[n - 1] * bin.inv(n);
        },
        n);
}

template<typename Tp> inline std::vector<Tp> fps_pow(std::vector<Tp> a, long long e, int n) {
    if (n <= 0) return {};
    if (e == 0) {
        std::vector<Tp> res(n);
        res[0] = 1;
        return res;
    }

    const int o = order(a);
    if (o < 0 || o > n / e || (o == n / e && n % e == 0)) return std::vector<Tp>(n);
    if (o != 0) a.erase(a.begin(), a.begin() + o);

    const Tp ia0 = a[0].inv();
    const Tp a0e = a[0].pow(e);
    const Tp me  = e;

    for (int i = 0; i < (int)a.size(); ++i) a[i] *= ia0;
    a = fps_log(a, n - o * e);
    for (int i = 0; i < (int)a.size(); ++i) a[i] *= me;
    a = fps_exp(a, n - o * e);
    for (int i = 0; i < (int)a.size(); ++i) a[i] *= a0e;

    a.insert(a.begin(), o * e, Tp(0));
    return a;
}
#line 7 "poly_basic.hpp"
#include <array>
#line 10 "poly_basic.hpp"

template<typename Tp> inline int degree(const std::vector<Tp> &a) {
    int n = (int)a.size() - 1;
    while (n >= 0 && a[n] == 0) --n;
    return n;
}

template<typename Tp> inline void shrink(std::vector<Tp> &a) { a.resize(degree(a) + 1); }

template<typename Tp> inline std::vector<Tp> taylor_shift(std::vector<Tp> a, Tp c) {
    const int n = a.size();
    auto &&bin  = Binomial<Tp>::get(n);
    for (int i = 0; i < n; ++i) a[i] *= bin.factorial(i);
    Tp cc = 1;
    std::vector<Tp> b(n);
    for (int i = 0; i < n; ++i) {
        b[i] = cc * bin.inv_factorial(i);
        cc *= c;
    }
    std::reverse(a.begin(), a.end());
    auto ab = convolution(a, b);
    ab.resize(n);
    std::reverse(ab.begin(), ab.end());
    for (int i = 0; i < n; ++i) ab[i] *= bin.inv_factorial(i);
    return ab;
}

// returns (quotient, remainder)
// O(deg(Q)deg(B))
template<typename Tp> inline std::array<std::vector<Tp>, 2>
euclidean_div_naive(const std::vector<Tp> &A, const std::vector<Tp> &B) {
    const int degA = degree(A);
    const int degB = degree(B);
    assert(degB >= 0);
    const int degQ = degA - degB;
    if (degQ < 0) return {std::vector<Tp>{Tp(0)}, A};
    std::vector<Tp> Q(degQ + 1), R = A;
    const auto inv = B[degB].inv();
    for (int i = degQ, n = degA; i >= 0; --i)
        if ((Q[i] = R[n--] * inv) != 0)
            for (int j = 0; j <= degB; ++j) R[i + j] -= Q[i] * B[j];
    R.resize(degB);
    return {Q, R};
}

// O(min(deg(Q)^2,deg(Q)deg(B)))
template<typename Tp> inline std::vector<Tp>
euclidean_div_quotient_naive(const std::vector<Tp> &A, const std::vector<Tp> &B) {
    const int degA = degree(A);
    const int degB = degree(B);
    assert(degB >= 0);
    const int degQ = degA - degB;
    if (degQ < 0) return {Tp(0)};
    const auto inv = B[degB].inv();
    std::vector<Tp> Q(degQ + 1);
    for (int i = 0; i <= degQ; ++i) {
        for (int j = 1; j <= std::min(i, degB); ++j) Q[degQ - i] += B[degB - j] * Q[degQ - i + j];
        Q[degQ - i] = (A[degA - i] - Q[degQ - i]) * inv;
    }
    return Q;
}

// returns (quotient, remainder)
template<typename Tp> inline std::array<std::vector<Tp>, 2>
euclidean_div(const std::vector<Tp> &A, const std::vector<Tp> &B) {
    const int degA = degree(A);
    const int degB = degree(B);
    assert(degB >= 0);
    // A = Q*B + R => A/B = Q + R/B in R((x^(-1)))
    const int degQ = degA - degB;
    if (degQ < 0) return {std::vector<Tp>{Tp(0)}, A};
    if (degQ < 60 || degB < 60) return euclidean_div_naive(A, B);

    auto Q = fps_div(std::vector(A.rend() - (degA + 1), A.rend()),
                     std::vector(B.rend() - (degB + 1), B.rend()), degQ + 1);
    std::reverse(Q.begin(), Q.end());

    // returns a mod (x^n-1)
    auto make_cyclic = [](const std::vector<Tp> &a, int n) {
        assert((n & (n - 1)) == 0);
        std::vector<Tp> b(n);
        for (int i = 0; i < (int)a.size(); ++i) b[i & (n - 1)] += a[i];
        return b;
    };

    const int len      = fft_len(std::max(degB, 1));
    const auto cyclicA = make_cyclic(A, len);
    auto cyclicB       = make_cyclic(B, len);
    auto cyclicQ       = make_cyclic(Q, len);

    fft(cyclicQ);
    fft(cyclicB);
    for (int i = 0; i < len; ++i) cyclicQ[i] *= cyclicB[i];
    inv_fft(cyclicQ);

    // R = A - QB mod (x^n-1) (n >= degB)
    std::vector<Tp> R(degB);
    for (int i = 0; i < degB; ++i) R[i] = cyclicA[i] - cyclicQ[i];
    return {Q, R};
}

template<typename Tp>
inline std::vector<Tp> euclidean_div_quotient(const std::vector<Tp> &A, const std::vector<Tp> &B) {
    const int degA = degree(A);
    const int degB = degree(B);
    assert(degB >= 0);
    // A = Q*B + R => A/B = Q + R/B in R((x^(-1)))
    const int degQ = degA - degB;
    if (degQ < 0) return {Tp(0)};
    if (std::min(degQ, degB) < 60) return euclidean_div_quotient_naive(A, B);

    auto Q = fps_div(std::vector(A.rend() - (degA + 1), A.rend()),
                     std::vector(B.rend() - (degB + 1), B.rend()), degQ + 1);
    std::reverse(Q.begin(), Q.end());
    return Q;
}
#line 8 "czt.hpp"

// returns F(a),F(ac),F(ac^2),...,F(ac^(n-1))
// Use        ij = binom(i,2)   + binom(-j,2) - binom(i-j,2)
// instead of ij = binom(i+j,2) - binom(i,2)  - binom(j,2)
template<typename Tp> inline std::vector<Tp> czt(std::vector<Tp> F, Tp c, int n, Tp a = 1) {
    if (n <= 0) return {};
    const int degF = degree(F);
    shrink(F);
    if (degF < 0) return std::vector<Tp>(n);
    if (degF == 0 || a == 0) return std::vector<Tp>(n, F[0]);
    if (a != 1) {
        // F(x) <- F(ax)
        Tp aa = 1;
        for (int i = 0; i <= degF; ++i) F[i] *= aa, aa *= a;
    }
    if (c == 0) {
        std::vector<Tp> res(n, F[0]);
        for (int i = 1; i <= degF; ++i) res[0] += F[i];
        return res;
    }
    std::vector<Tp> H(std::max(degF + 1, n - 1));
    Tp cc = H[0] = 1;
    for (int i = 1; i < (int)H.size(); ++i) H[i] = H[i - 1] * (cc *= c);
    std::vector<Tp> G(degF + n); // G[i+degF]=c^(-binom(i,2))
    auto GG     = G.begin() + degF;
    const Tp ic = c.inv();
    cc = GG[0] = 1;
    for (int i = 1; i < n; ++i) GG[i] = GG[i - 1] * cc, cc *= ic;
    cc = 1;
    for (int i = -1; i >= -degF; --i) GG[i] = GG[i + 1] * (cc *= ic);
    // F[i] <- c^(binom(i+1,2))*F[i]
    for (int i = 0; i <= degF; ++i) F[i] *= H[i];
    F = middle_product(G, F);
    // F[i] <- c^(binom(i,2))*F[i]
    for (int i = 1; i < n; ++i) F[i] *= H[i - 1];
    return F;
}

// returns f s.t. f(aq^i)=F[i]
// aq^i != aq^j for all i != j
// see: https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric
// noshi91. 標本点が等比数列を成す場合に補間多項式を計算するアルゴリズム.
template<typename Tp> inline std::vector<Tp> inv_czt(const std::vector<Tp> &F, Tp q, Tp a = 1) {
    if (F.empty()) return {};
    if (a == 0) return {F[0]};
    const int n = F.size();
    std::vector<Tp> Q(n), S(n), M(n), D(n);
    Tp qq = 1;
    // Q[i]=q^i
    for (int i = 0; i < n; ++i) Q[i] = qq, qq *= q;
    // S[i]=prod_(i=1..i)(1-q^i)
    S[0] = 1;
    for (int i = 1; i < n; ++i) S[i] = S[i - 1] * (1 - Q[i]);
    const auto Sn   = S[n - 1] * (1 - qq);
    const auto invS = batch_inv(S);
    qq              = 1;
    // M[n-i]=qbinom(n,i)*q^(binom(i,2))*(-1)^i
    for (int i = 1; i < n; ++i) M[n - i] = Sn * invS[i] * invS[n - i] * (qq *= -Q[i - 1]);
    M[0] = qq * -Q[n - 1]; // in case of q^n=1
    // D[i]=S[i]*S[n-i-1]*q^(binom(i,2)+i(n-i-1))*(-1)^i
    D[0] = 1;
    for (int i = 0; i < n - 1; ++i) D[i + 1] = D[i] * -Q[n - i - 2];
    for (int i = 0; i < n; ++i) D[i] *= S[i] * S[n - i - 1];
    // D[i] <- -F[i]/D[i]
    D = batch_inv(D);
    for (int i = 0; i < n; ++i) D[i] *= -F[i];
    auto res = convolution(M, czt(D, q.inv(), n, q.inv()));
    res.resize(n);
    if (a != 1) {
        const auto ia = a.inv();
        Tp aa         = 1;
        for (int i = 0; i < n; ++i) res[i] *= aa, aa *= ia;
    }
    return res;
}

Chirp Z Transform (czt.hpp)