Chirp-Z Transform / Inverse Chirp-Z Transform

· · 算法·理论

Chirp-Z Transform / Inverse Chirp-Z Transform

读者需要具有一定的多项式基础和数学基础。

Chirp-Z Transform

给定 n-1 次多项式 F(x) ，求它在每个 i\in[0,n-1],p\cdot q^i 上的点值，要求 \mathbf O(n\log n)。

Chirp-Z Transform（CZT）解决的是上面的问题，概括出来就是等比数列多项式多点求值。

我们先尝试构造一个卷积的形式：

y_i=\sum_{j=0}^{n-1}f_j\cdot x^j=\sum_{j=0}^{n-1}f_j(p\cdot q^i)^j=\sum_{j=0}^{n-1}(f_j\cdot p^j)q^{i\cdot j}

解释：根据点值的定义，将 x_i=p\cdot q^i 代入，然后变形将指数 j 给到括号内，提出 q 相关的项。

y_i=\sum_{j=0}^{n-1}(f_j\cdot p^j)q^{\binom{i+j}2-\binom i2-\binom j2}

解释：有恒等式 i\cdot j=\binom{i+j}2-\binom i2-\binom j2。

直接按组合数的定义展开即证。

\begin{aligned} &\binom{i+j}2-\binom i2-\binom j2\\ &=\frac12(i+j)(i+j-1)-\frac12i(i-1)-\frac12j(j-1)\\ &=\frac12(i^2+i\cdot j-i+i\cdot j+j^2-j-i^2+i-j^2+j)\\ &=\frac12(2\cdot i\cdot j)\\ &=i\cdot j \end{aligned}
y_i=q^{-\binom i2}\sum_{j=0}^{n-1}(f_j\cdot p^j)q^{\binom{i+j}2-\binom j2}=q^{-\binom i2}\sum_{j=0}^{n-1}(f_j\cdot p^j\cdot q^{-\binom j2})q^{\binom{i+j}2}

解释：将与求和索引 j 无关的项提出到求和符号外，然后展开乘方运算，将仅与 j 有关的量放到一起，构造卷积形式。

y_i=q^{-\binom i2}\sum_{j=0}^{n-1}u_j\cdot v_{i+j}

解释：换元 u_t=f_t\cdot p^t\cdot q^{-\binom t2},v_t=q^{\binom t2}。此时后面的求和部分 \sum u_j\cdot v_{i+j} 是一个减法卷积形式，可以构造 U(x)=\sum u_j\cdot x^j,V(x)=\sum v_j\cdot x^j，计算 U(x) 和 V(x) 的减法卷积即可得到 \sum u_j\cdot v_{i+j} 部分，最后乘上前面的系数 q^{-\binom i2} 就完成了。

那如何计算 U(x) 和 V(x) 的减法卷积呢？经典的方法是将其中一个系数翻转，然后计算二者的加法卷积（多项式乘法）。使用快速数论变换实现多项式乘法可以做到 \mathbf O(n\log n) 的时间复杂度。

Inverse Chirp-Z Transform

给定 n-1 次多项式 F(x) 在每个 i\in[0,n-1],p\cdot q^i 上的点值，求它，要求 \mathbf O(n\log n)。

Inverse Chirp-Z Transform（ICZT）解决的是上面的问题，概括出来就是等比数列多项式多点插值。

我们先尝试简化问题。设 G(x)=\sum\limits_{j=0}^{n-1}g_j\cdot x^j=F(p\cdot x)=\sum\limits_{j=0}^{n-1}f_j(p^j\cdot x^j)。如果我们能求出 G(x) 的系数 g_j，又因为 g_j=f_j\cdot p^j，那我们就能直接得到 F(x) 的系数 f_j=g_j\cdot p^{-j}。此时问题简化为“给定 n-1 次多项式 G(x) 在每个 i\in[0,n-1],q^i 上的点值，求它，要求 \mathbf O(n\log n)”。

设 G(x) 在牛顿基下的形式为 G(x)=\sum\limits_{j=0}^{n-1}g'_j\prod\limits_{k=0}^{j-1}(x-q^k)。

y_i=g(q^i)=\sum_{j=0}^{n-1}g'_j\prod\limits_{k=0}^{j-1}(q^i-q^k)=\sum_{j=0}^ig'_j\prod_{k=0}^{j-1}(q^i-q^k)

解释：根据 G(x) 的牛顿基表示形式展开，然后根据点值的定义带入。后面我们把第一个求和式中 j 的上界从 n-1 缩小为 i，是因为若 j>i，则第二个求积式中 k 可以取到 i，导致乘积结果为 0，对 y_i 无影响。

y_i=\sum_{j=0}^ig'_j\cdot q^{\binom j2}\dfrac{(i)_q!}{(i-j)_q!}

解释：根据变形 q-阶乘恒等式 \prod\limits_{k=0}^{j-1}(q^i-q^k)=q^{\binom j2}\dfrac{(i)_q!}{(i-j)_q!} 可得。

请不要被 q-阶乘这个名字吓住。这只是个定义：(n)_q!=\prod\limits_{i=1}^n(q^i-1)。熟悉 q-analog 的读者应该会发现，这不是标准的 q-阶乘形式 [n]_q!。这里是我们有意为之，定义此变形 q-阶乘，舍去了分母的部分，并且将分子的部分乘上了 -1。

\prod_{k=0}^{j-1}(q^i-q^k)=\prod_{k=0}^{j-1}q^k(q^{i-k}-1)=q^{\sum\limits_{k=0}^{j-1}k}\prod_{k=0}^{j-1}(q^{i-k}-1)

解释：提出 q^k 一项，然后通过乘方运算规律将其完全扯到求积式外面。

\prod_{k=0}^{j-1}(q^i-q^k)=q^{\binom j2}\prod_{k=0}^{j-1}(q^{i-k}-1)

解释：求积式外面的系数的指数是等差数列求和，其结果为 \frac{j(j-1)}2，这等同于 \binom j2。观察目前这个式子与目标式 q^{\binom j2}\dfrac{(i)_q!}{(i-j)_q!} 的差异，接下来我们只需要证明 \prod\limits_{k=0}^{j-1}(q^{i-k}-1)=\dfrac{(i)_q!}{(i-j)_q!} 即可。

\dfrac{(i)_q!}{(i-j)_q!}=\prod_{l=i-j+1}^i(q^l-1)

解释：根据定义有 (i)_q!=\prod\limits_{l=1}^i(q^l-1) 以及 (i-j)_q!=\prod\limits_{l=1}^{i-j}(q^l-1)，直接相除即可。

\dfrac{(i)_q!}{(i-j)_q!}=\prod_{k=0}^{j-1}(q^{i-k}-1)

解释：换元 k=i-l，也就有 l=i-k。原先 l=i-j+1 的下界对应现在 i-k=i-j+1 也就是 k=j-1 的上界；原先 l=i 的上界对应现在 i-k=i 也就是 k=0 的下界，证毕。

如果想了解更多关于 q-analog 的内容可以阅读 WorldMachine 的专栏 - 分拆数与 q-analog 小记。

\dfrac{y_i}{(i)_q!}=\sum_{j=0}^i\left(g'_j\cdot q^{\binom j2}\right)\dfrac1{(i-j)_q!}

解释：移项整理，尝试构造卷积形式。

a_i=\sum_{j=0}^ib_j\cdot c_{i-j}

解释：换元 a_t=\dfrac{y_t}{(t)_q!},b_t=g'_t\cdot q^{\binom t2},c_t=\dfrac1{(t)_q!}。此时后面的求和部分 \sum b_j\cdot c_{i-j} 是一个加法卷积形式，可以构造 B(x)=\sum b_j\cdot x^j,C(x)=\sum c_j\cdot x^j。注意你现在要求的是 g'_j 一项，也就是要求 B(x)，那就需要多项式除法（多项式乘法逆元）求出 B(x)，然后 g'_j=b_j\cdot q^{-\binom j2} 即可。

接下来我们需要把牛顿基系数 g' 转换回常规系数 g。

G(x)=\sum_{j=0}^{n-1}g'_j\prod_{k=0}^{j-1}(x-q^k)=\sum_{j=0}^{n-1}g'_j\sum_{k=0}^j\binom jk_q(-1)^{j-k}q^{\binom{j-k}2}x^k

解释：牛顿基式子，然后 q-二项式定理把后面的求积式化掉。

没想到吧 q-analog 还在发力。q-二项式定理的内容为：

\prod_{k=0}^{j-1}(x-q^k)=\sum_{k=0}^j\binom jk_q (-1)^{j-k}q^{\binom{j-k}2}x^k

首先定义 q-二项式系数 \binom nk_q=\frac{(n)_q!}{(k)_q!(n-k)_q!}。注意这里的 (n)_q! 是我们之前定义的变形 q-阶乘 \prod\limits_{i=1}^n(q^i-1)。我们需要先证明一个递推关系 \binom nk_q=\binom{n-1}{k-1}_q+q^k\binom{n-1}k_q，它又名 q-帕斯卡恒等式。

\binom{n-1}{k-1}_q+q^k\binom{n-1}k_q=\frac{(n-1)_q!}{(k-1)_q!(n-k)_q!}+q^k\frac{(n-1)_q!}{(k)_q!(n-1-k)_q!}

解释：按照 q-二项式系数的定义展开。

\binom{n-1}{k-1}_q+q^k\binom{n-1}k_q=\frac{(n-1)_q!}{(k)_q!(n-k)_q!}\left((q^k-1)+q^k(q^{n-k}-1)\right)

解释：通分。

\binom{n-1}{k-1}_q+q^k\binom{n-1}k_q=\frac{(n-1)_q!}{(k)_q!(n-k)_q!}(q^n-1)

解释：展开分式右侧的括号内容。

\binom{n-1}{k-1}_q+q^k\binom{n-1}k_q=\frac{(n)_q!}{(k)_q!(n-k)_q!}=\binom nk_q

解释：根据定义有 (n)_q!=(n-1)_q!(q^n-1)，直接把括号内容乘到分子上合并。然后得到的式子根据 q-二项式系数的定义变形，证毕。

接下来我们使用数学归纳法证明 q-二项式定理，先是一些记号：

  • 左 L_n(x)=\prod\limits_{k=0}^{n-1}(x-q^k)，特别地 L_n(0) 是空积，也就是说 L_n(0)=1。
  • 右 R_n(x)=\sum\limits_{k=0}^nc_{n,k}x^k，这里我们换元 c_{n,k}=\binom nk_q(-1)^{n-k}q^{\binom {n-k}2} 方便书写。

下面开始证明。当 n=0 时，L_n(x)=1，R_n(x)=\binom 00_q(-1)^0q^0x^0=1，结论成立。假设对 n-1 成立，我们知道 L_n(x)=L_{n-1}(x)(x-q^{n-1})，考虑这个等式左右两边每一项的关系：

\forall k,[x^k]L_n(x)=[x^{k-1}]L_{n-1}(x)-q^{n-1}[x^k]L_{n-1}(x)

解释：[x^k]L_n(x) 表示提取多项式 L_n(x) 的 x^k 项系数。根据前面提到的 L_n(x) 与 L_{n-1}(x) 之间的关系，我们知道 L_n(x) 的 x^k 项是由 L_{n-1} 的 x^{k-1} 项乘上 x 加上 L_{n-1} 的 x^k 项乘上 -q^{n-1} 得到。

\forall k,[x^k]L_n(x)=c_{n-1,k-1}-q^{n-1}c_{n-1,k}

解释：根据归纳假设，L_{n-1}(x)=R_{n-1}(x)，而 R_{n-1}(x) 是一个求和式，它的 x^k 项就是 c_{n-1,k}，同理它的 x^{k-1} 项就是 c_{n-1,k-1}。

\forall k,[x^k]L_n(x)=\binom{n-1}{k-1}_q(-1)^{n-k}q^{\binom{n-k}2}-q^{n-1}\binom{n-1}k_q(-1)^{n-1-k}q^{\binom{n-1-k}2}

解释：按照 c_{n,k} 的定义直接展开。

\forall k,[x^k]L_n(x)=(-1)^{n-k}q^{\binom{n-k}2}\left(\binom{n-1}{k-1}_q-q^{n-1}\binom{n-1}k_q(-1)^{-1}q^{\binom{n-1-k}2-\binom{n-k}2}\right)

解释：提取公因式。

\forall k,[x^k]L_n(x)=(-1)^{n-k}q^{\binom{n-k}2}\left(\binom{n-1}{k-1}_q+q^k\binom{n-1}k_q\right)

解释：首先将 (-1)^{-1} 直接乘掉。然后右边 q 的指数是 \binom{n-1-k}2-\binom{n-k}2=\frac{(n-k-1)(n-k-2)}2-\frac{(n-k)(n-k-1)}2=\frac{n-k-1}2(n-k-2-(n-k))=\frac{n-k-1}2(-2)=-n+k+1，前面还乘了一项 q^{n-1}，那么 q 的总指数就是 -n+k+1+n-1=k。

\forall k,[x^k]L_n(x)=(-1)^{n-k}q^{\binom{n-k}2}\binom nk_q=c_{n,k}=[x^k]R_n(x)

解释：第一步是 q-帕斯卡恒等式，然后根据 c_{n,k} 的定义直接化简。然后 c_{n,k} 又是 R_n(x) 的 x^k 项系数，所以我们证明了对于每个 k，L_n(x) 的 x^k 项系数和 R_n(x) 的 x^k 项系数相同，那么就有 L_n(x)=R_n(x)，q-二项式定理证毕。

交换求和顺序，对比系数就有 g_j=\sum\limits_{i=j}^{n-1}g'_i\binom ij_q(-1)^{i-j}q^{\binom{i-j}2}。

g_j=\sum_{i=j}^{n-1}g'_i\dfrac{(i)_q!}{(j)_q!(i-j)_q!}(-1)^{i-j}q^{\binom{i-j}2}

解释：根据 q-组合数的定义展开。

g_j(j)_q!=\sum_{i=j}^{n-1}(g'_i(i)_q!)\dfrac{(-1)^{i-j}q^{\binom{i-j}2}}{(i-j)_q!}

解释：移项整理，构造卷积形式。

g_j=\dfrac1{(j)_q!}\sum_{i=j}^{n-1}u_i\cdot v_{i-j}

解释：换元 u_t=g'_t(t)_q!,v_t=\dfrac{(-1)^tq^{\binom t2}}{(t)_q!}。此时后面的求和部分 \sum u_i\cdot v_{i-j} 是一个加法卷积形式，可以构造 U(x)=\sum u_j\cdot x^j,V(x)=\sum v_j\cdot x^j。计算 U(x) 和 V(x) 的加法卷积即可得到 \sum u_j\cdot v_{i-j} 部分，最后乘上前面的系数 \dfrac1{(j)_q!} 就完成了。然后 f_j=g_j\cdot p^{-j}，得到 f。

以上多项式乘法逆元可以牛顿迭代，加法卷积可以快速数论变换，可以做到 \mathbf O(n\log n) 的时间复杂度。