如何在任意代数结构上做多项式乘法

前两天水群的时候看到有人说可以模 \(2^{64}\) 做多项式乘法且不用 MTT 让我属实震撼，然后去翻了翻文献找到了做法。

做法来自 Cantor 在 1991 年发的论文，这个算法也叫作 Cantor's Algorithm。

全名是 David Geoffrey Cantor，不是那个集合论的 Cantor。

我们解决了什么

根据原文献的说法，我们可以做到 "multiplying polynomials with coefficients from an arbitrary, not necessarily commutative, not necessarily associative, algebra \(\mathcal A\)"。

这个「任意」挺让人摸不着头脑的，加上后面的「非交换，非结合」就有种营销号的感觉。然后后面文献详细做出了解释：

Specifically the algebra \(\mathcal A\) must be an Abelian group under "\(+\)" and have a binary operation "\(\cdot\)" satisfying the distributive law

\[(u+v)\cdot(x+y)=u\cdot x+u\cdot y+v\cdot x+v\cdot y \]

换句话说，就是要求 \(\mathcal A\) 是个不用满足乘法结合律的环，好像一般也叫作「非结合环」。

一个例子是我们现在可以对四元数/矩阵为系数的多项式做乘法。

关于复杂度：做两个度为 \(n\) 的多项式乘法，需要做 \(\mathcal O(n\log n)\) 次 \(\mathcal A\) 中的乘法和 \(\mathcal O(n\log n\log\log n)\) 次 \(\mathcal A\) 中的加法，对某些特殊的代数结构可以做特化让复杂度变成 \(\mathcal O(n\log n)\)。

算法原理

考虑一下我们在做多项式乘法的时候用到了代数结构的什么性质：

• 假设我们要做长为 \(n\) 的 IDFT，则 \(n\) 必须可逆（因为在 IDFT 的时候用到了除法）
• 我们需要 \(n\) 阶单位根 \(\omega_n\)，且 \(\sum_{j=0}^{n-1}\omega_n^{ij}=[i=0]n\)。

我们分三部分解决问题。

Part 1：如何规避除法

下文中的代数结构中元素和整数的乘法的含义是若干个代数结构中元素的和。

假设现在要做 \(n=s^r\) 阶的 DFT，其中不要求 \(s\) 是素数，而且代数结构内有 \(ns\) 阶单位根 \(\omega_{ns}\)。

设要做多项式 \(A(x)=\sum_{i=0}^{n-1}a_ix^i\) 和 \(B(x)=\sum_{i=0}^{n-1} b_ix^i\) 的乘法，结果为 \(C(x)=\sum_{i=0}^{2n-2} c_ix^i\)。设 \(\{A_i\}\) 和 \(\{B_i\}\) 是它们 DFT 后的数组，然后令 \(D_i=A_iB_i\)，那么考虑对 \(\{D_i\}\) 做 IDFT 后还没做除法前的结果为 \(\{d_i\}\)，那么有

\[d_h=n(\sum_{i\equiv h\pmod n}c_i)=n(c_h+c_{h+n}) \]

对 \(\{a_i\}\) 做变换变成 \(\{\hat a_i=a_i\omega_{ns}^i\}\)，对 \(\{b_i\}\) 做同样的变换变成 \(\{\hat b_i\}\)，对它们做 DFT 变成 \(\{\hat A_i\}\) 和 \(\{\hat B_i\}\)，然后令 \(E_i=\hat A_i\hat B_i\)，设 \(\{e_i\omega_{ns}^i\}\) 为对 \(\{E_i\}\) 做 DFT 还没做除法的结果，那么可以得到：

\[e_h\omega_{ns}^h=n(\sum_{i\equiv h\pmod n}c_i\omega_{ns}^i)=n(c_h\omega_{ns}^h+c_{h+n}\omega_{ns}^{h+n})\\ \implies e_h=n(c_h+c_{h+n}\omega_s) \]

于是得到

\[\left\{ \begin{aligned} &n(1-\omega_s)c_h=e_h-d_h\omega_s\\ &n(1-\omega_s)c_{n+h}=d_h-e_h \end{aligned} \right. \]

设

\[\tau_n=\prod_{1\leq i<n,\gcd(i,n)=1}(1-\omega_n^i) \]

那么根据某个结论，当 \(n\) 为素数 \(p\) 的幂的时候 \(\tau_n=p\)，否则 \(\tau_n=1\)，它一定是个整数。这蕴含了

\[\frac{\tau_s}{1-\omega_s}=\prod_{2\leq i<s,\gcd(i,s)=1}(1-\omega_s^i) \]

于是

\[\left\{ \begin{aligned} &n\tau_sc_i=(e_i-d_i\omega_s)\prod_{2\leq i<s,\gcd(i,s)=1}(1-\omega_s^i)\\ &n\tau_sc_{i+n}=(d_i-e_i)\prod_{2\leq i<s,\gcd(i,s)=1}(1-\omega_s^i) \end{aligned} \right. \]

我们随便选两个互素的 \(s\) 做这个东西，就可以得到 \(M_1c_i\) 和 \(M_2c_i\) 的值，然后注意到 \(M_1\) 和 \(M_2\) 是互素的，于是可以对它们做扩展欧几里得得到 \(k_1M_1+k_2M_2=1\)，然后就可以在不做除法的前提下得到 \(c_i\) 的值。

Lemma 1. 如果存在单位根，则存在算法，只需要进行

• \(\mathcal O(rn^2s\mu_s)\) 次加法
• \(\mathcal O(n)\) 次乘法

即可计算出所有 \(\tau_snc_i\)。

其中 \(\mu_s\) 是 \(\Phi_s\) 中的所有系数的绝对值之和。

Part 2：如何规避单位根

做法非常简单，既然没有单位根我们就自己造！

先考虑一个不成熟的想法：类似扩域做，我们暴力引入 \(n\) 阶单位根 \(\omega_n\) 并带着它做运算，也就是说我们暴力把代数结构从 \(\mathcal A\) 变成了 \(\mathcal I_n=\mathcal A[\omega_n]\)。这样做的复杂度显然不低于 \(\mathcal O(n^2)\)。

然后考虑优化，先考虑一下 \(\omega_n\) 的实质是什么。

这里有个结论：\(\mathcal I_n\cong\mathcal A[x]/\lang\Phi_n(x)\rang\)，其中 \(\Phi_n(x)\) 是 \(n\) 阶分圆多项式，定义为

\[\Phi_n(x)=\prod_{0\leq i<n,\gcd(i,n)=1}(x-\omega_n^i) \]

换句话说，两个带着 \(\omega_n\) 的东西的运算相当于对应的多项式运算最后对分圆多项式 \(\Phi_n(x)\) 取模。

Example. \((1+\omega_3)^2=\omega_3\)

Proof.

上式显然是 \(1+2\omega_3+\omega_3^2\)。而注意到 \(\omega_3^2=-1-\omega_3\)，于是就是 \(\omega_3\)。

考虑成多项式乘法，\(\Phi_3(x)=1+x+x^2\)，于是相当于 \(1+2x+x^2\bmod (1+x+x^2)\)，这显然是 \(x\)。

可以证明，分圆多项式都是整系数素多项式，于是它一定是 \(\omega_n\) 的最小多项式，于是上述结论是不难证明的。

分圆多项式有性质：\(\Phi_n(x)|(x^n-1)\)，于是对 \(\Phi_n(x)\) 取模可以考虑成先对 \(x^n-1\) 取模，然后再 reduce 到 \(\Phi_n(x)\)。

我们后面会说明为什么不直接对 \(x^n-1\) 取模而是对 \(\Phi_n(x)\) 这种阴间东西取模。

接下来有个引理。

Lemma 2. 假设 \(n=s^r,N=s^u,u\geq r\)，则对一个系数在 \(\mathcal I_N\) 中的多项式做长为 \(n\) 的 DFT 的复杂度不超过 \(\mathcal O(snNr)\)。

Proof. 考虑一下做 DFT 的过程，先考虑 \(s=2\)，我们每次做的操作是

\[F(\omega_n^k)=F_0(\omega_{n/2}^k)+\omega_n^kF_1(\omega_{n/2}^k)\\ F(\omega_n^{k+n/2})=F_0(\omega_{n/2}^k)-\omega_n^kF_1(\omega_{n/2}^k) \]

\(s\neq 2\) 的操作是类似的，于是做的操作只有加法和乘以 \(\omega_N^k\)，而后者由于是对 \(x^N-1\) 取模，所以可以考虑成循环移位，于是两者都可以 \(\mathcal O(N)\) 做。总复杂度是 \(T(n)=sT(\frac ns)+\mathcal O(snN)=\mathcal O(snNr)\)。

最后答案对 \(\Phi_n(x)\) 取模的过程不会影响这个复杂度。\(\blacksquare\)

接下来是重头戏：

Lemma 3. 如果可以进行除法，存在算法，做两个度不超过 \(n\) 的多项式 \(A(x)\) 和 \(B(x)\) 的乘法，需要不超过 \(\mathcal O(n\log n)\) 次对应代数中的乘法和不超过 \(\mathcal O(n\log n\log\log n)\) 次对应代数中的加法。

Proof.

设 \(m=s^r\) 满足 \(\varphi(m)>\deg A+\deg B\)。我们不考虑做多项式乘法 \(A(x)B(x)\)，而是直接考虑成计算 \(A(\omega_m)B(\omega_m)\)。

注意到答案是对 \(\Phi_m(x)\) 取模，而 \(\Phi_m(x)\) 是素多项式，于是不会出现两个不相等的多项式等价。而又由于 \(\deg \Phi_m(x)=\varphi(m)>\deg A+\deg B\)，所以这个转化是完全合法的！只需要最后把每一项从 \(\omega_m\) 转化成 \(x\) 即可得到答案。

换言之，我们现在在对 \(\mathcal I_m[x]/\lang x-\omega_m\rang\) 上的元素做乘法，也就是说 \(x\) 和 \(\omega_m\) 是等价的元素。

设 \(p=s^u,q=s^v\) 且 \(u+v=r,v+1\leq u\leq v+2\)，然后我们把多项式写成这样的形式。

\[\begin{aligned} A(x)&=\sum_{j=0}^{q-1}(\sum_{i=0}^{\varphi(p)-1}a_{iq+j}x^{iq})x^j\\ &=\sum_{j=0}^{q-1}(\sum_{i=0}^{\varphi(p)-1}a_{iq+j}\omega_m^{iq})x^j\\ &=\sum_{j=0}^{q-1}(\sum_{i=0}^{\varphi(p)-1}a_{iq+j}\omega_p^{i})x^j\\ \end{aligned} \]

现在是一个神仙转化：我们把 \(A(x)\) 看成 \(\mathcal I_p[x]\) 中度为 \(q\) 的元素再进行 DFT，也就是说，我们把中间那一坨东西看成系数然后进行长为 \(q\) 的 DFT。也就是说，我们通过一些转换，凭空创造出了单位根！

考虑一下这个东西的复杂度。由于 \(q\leq p\)，所以可以使用 Lemma 进行分析，进行 DFT 的复杂度不超过 \(\mathcal O(spqr)=\mathcal O(snr)\)。而我们使用充分小的 \(s\) 进行 DFT，于是这个复杂度是 \(\mathcal O(n\log n)\)，其中只进行了加法。

做完 DFT 之后我们需要计算两个 \(\mathcal I_p\) 中元素的乘法，这可以递归调用上述过程。

加法的复杂度分析：\(T(n)=\sqrt nT(\sqrt n)+\mathcal O(n\log n)=\mathcal O(n\log n\log\log n)\)。

乘法的复杂度分析：\(T(n)=\sqrt nT(n)+\mathcal O(1)=\mathcal O(n\log n)\)。

\(\blacksquare\)

具体分析可以参见原论文。

Part 3：如何同时规避单位根和除法

发现 Part 2 直接为我们凭空创造出了单位根，直接拿着在 Part 1 当中用即可。

一般会带入 \(s=2\) 和 \(s=3\) 进行计算，于是所有和 \(s\) 有关的东西都可以看成常数，于是复杂度可以直接 merge 起来。

Theorem. 对任意代数结构 \(\mathcal A\)，存在算法，计算两个度为 \(n\) 且系数为 \(\mathcal A\) 中元素的多项式只需要 \(\mathcal O(n\log n\log \log n)\) 次 \(\mathcal A\) 上的加法和 \(\mathcal O(n\log n)\) 次 \(\mathcal A\) 上的乘法。

参考文献

[1] Fast convolution for 64-bit integers

[2] D. G. Cantor and E. Kaltofen: On Fast Multiplication of Polynomials over Arbitrary Algebra, 1991

[3] 彭雨翔: Introduction to Polynomials

注：[3] 可以在 U 群找到。