Howto Multiply Great Numbers

发表于 2026/01/08 更新于 2026/01/22

作者 Li 6in9Muyou

16 分钟阅读

超级大数的乘法

我的问题是把两个特别大的数字相乘比如$1234\times5678$ 小学乘法需要逐个数字互相相乘用时是$O(n^2)$的这是显然的基本操作是较小数字的乘法比如上面例子中每个数字互相乘就是十以内乘法是较小数字的乘法

一个比较快的算法大体思路如下

第一先变换把两个数字的每一位看成一个数组的元素例如num.split('') 然后把这些数组做某个变换变成另一个数组这个数组里面的元素都是较小的数字肯定比两个乘数本身小

第二把这两个数组里的元素逐个相乘得到一个结果数组这个步骤是简单直白的

第三逆变换结果数组得到最终的结果

首先必须把数字写成多项式 $A(x)=1+2x+3x^2+4x^3$ 有多少位数就有多少个多项式系数两个多项式相乘乘积结果的系数个数是原先两个乘数的系数个数的和比如3位数乘5位数结果最多有8位数

先来看变换过程变换过程就是取足够多的数字放到每个乘数对应的多项式里面取求值足够多的意思就是有N个系数的多项式至少需要取N+1个数字那么逆变换过程就是解方程下面演示有2个系数的多项式相乘也就是$12\times45=540$

\[A(x) = 0x^2+1x + 2\\ B(x) = 0x^2+4x + 5\\\]

因为乘积有3位数变换时我们任取4个数字比如11、22、33、44 变换并逐对相乘之后得到

\[C(11) = A(11) \times B(11) = 13 \times 49 = 637\\ C(22) = A(22) \times B(22) = 24 \times 93 = 2232\\ C(33) = A(33) \times B(33) = 35 \times 137 = 4795\\ C(44) = A(44) \times B(44) = 46 \times 181 = 8326\]

逆变换过程就是解方程求解过程从略求得$C(x)=4x^2 + 13x + 10$

这个是不符合十进制的位数的需要进位进位后得到540

显然如果这几个数字是任意取的话是没有任何意义的变换过程和最后解方程过程是非常费力的解方程过程的费力程度先略过了在变换过程中用秦九韶算法对有2N+1个系数的多项式求值需要2N次乘法求值2N+2次就会需要N的平方次乘法跟小学乘法一样这样搞是没有意义的

所以需要取一些非常特殊的数字来简化这个多项式求值的过程怎么选呢选什么样的数字呢？

注意到$x^n$跟系数相乘的这个乘法是不可能省略的所以考虑在计算幂次的时候省略我们知道平方是有对称性的例如$(-2)^2=2^2=4$ 那么我们可以取一套相反数如在下例中 $123\times 456=56088$ 结果是5位数需要5个系数的多项式算出这些系数需要6个数不妨取77、88、99 和他们的相反数注意到如果按照$C(x)$的形式分解 $E, O$都只需要代入3次可以大致节省一半的乘法次数

\[A(x) = 1x^2+2x+3 = (1x^2 + 3) + x \cdot (2) = A_E(x^2) + x \cdot A_O(x^2)\\ B(x) = 4x^2+5x+6 = (4x^2 + 6) + x \cdot (5) = B_E(x^2) + x \cdot B_O(x^2)\\ C(x)=c_4x^4+c_3x^3+c_2x^2+c_1x+c_0=E(x^2)+x\cdot O(x^2)\]

也就是说可以利用平方的对称性利用方法就是把多项式拆成奇数次和偶数次两个部分每次拆分之后都可以节省一半的乘法次数如果能不断拆分下去还可以节省更多乘法次数但是拆分一次之后再往下拆分就不具有对称性了例如

\[D(x) = d_4x^4+d_3x^3+d_2x^2+d_1x+c_0 \\ = (d_4u^2 + d_2u + d_0) + x \cdot (d_3u + d_1)\\ = [(d_4w + d_0) + u \cdot (d_2)]+ x \cdot [(d1)+u \cdot (d_3)]\]

在上例中三组相反数就不能为$w$项节省乘法次数了用$w$计算跟用$u^2$计算一样都是需要代入3次

对于多项式 $D(x) = d_4x^4 + d_3x^3 + d_2x^2 + d_1x + d_0$，我们可以按照次数的奇偶性将其项进行归类。这种拆分是 快速傅里叶变换 (FFT) 分治步骤的核心逻辑：

1. 拆分过程

我们观察各项的幂次：

偶次项：$d_4x^4, d_2x^2, d_0$
奇次项：$d_3x^3, d_1x$

我们将 $D(x)$ 写成这两个部分的和：

2. 引入 $E(x^2)$ 和 $O(x^2)$

为了让分治递归能够进行，我们需要将变量统一化，即让括号内的多项式都变成关于 $x^2$ 的函数：

偶数部分 $E(x^2)$ (Even part)：
\[E(x^2) = d_4(x^2)^2 + d_2(x^2)^1 + d_0(x^2)^0\]
令 $z = x^2$，则 $E(z) = d_4z^2 + d_2z + d_0$。
奇数部分 $O(x^2)$ (Odd part)：
从奇次项中提取一个 $x$：
\[d_3x^3 + d_1x = x(d_3x^2 + d_1)\]
令 $O(x^2) = d_3(x^2)^1 + d_1(x^2)^0$。
同样令 $z = x^2$，则 $O(z) = d_3z + d_1$。

3. 最终分解形式

将上述部分代回原式，得到：

\[D(x) = E(x^2) + x \cdot O(x^2)\]

4. 为什么在大数乘法里这很重要？

当我们代入一对相反数 $x$ 和 $-x$ 时，由于 $(x)^2 = (-x)^2 = z$，我们只需要计算一次 $E(z)$ 和 $O(z)$ 的值，就可以直接通过一次加法和一次减法算出两个点值：

\[\begin{cases} D(x) = E(z) + x \cdot O(z) \\ D(-x) = E(z) - x \cdot O(z) \end{cases}\]

进阶工况思考：

对于规模极大的输入（如千万位的大数），$E(x^2)$ 和 $O(x^2)$ 本身又可以作为新的多项式继续向下进行奇偶拆分，直到多项式只剩下一个常数项。这种递归过程将原本 $O(N^2)$ 的复杂度降低到了 $O(N \log N)$。

在硬件实现层面，这种分解允许我们使用 SIMD (单指令多数据流) 指令集（如 AVX-512）同时处理 $x$ 和 $-x$ 的计算，极大地提高了吞吐量。

接下来，你想看一看如何对 $E(z)$ 进行进一步的递归拆分，还是想通过具体的数值（比如代入 $x=77$）来验证这个公式的等价性？

单位根是指满足等式 $x^n = 1$ 的所有复数（或在模空间中的大整数）。它们具有以下四个支撑起整个快速变换算法的核心性质：

1. 周期性（Periodicity）

单位根是循环的。如果你在圆周上一直转，最终会回到起点。

数学表示：$\omega_n^{k+n} = \omega_n^k$
工业意义：在处理数量特别大的输入时，我们不需要存储无限个幂次。对于 $N=1048576$ 的变换，所有的旋转因子都落在 $[0, N-1]$ 这个有限的集合里。这使得硬件预取（Pre-fetching）旋转因子表变得极其高效。

2. 折半性（Halving Property） —— 分治算法的基石

这是 NTT 能够通过奇偶拆分（EO 拆分）不断降维的根本原因。

数学表示：$(\omega_n^k)^2 = \omega_{n/2}^k$
深度拆解：
一个 $N$ 次单位根的平方，恰好等于 $N/2$ 次单位根。
演示：
在 $N=8$ 的求值中，代入 8 个点。当你把这 8 个点平方后（即处理偶数项和奇数项时），由于这个性质，它们变成了 4 个重复的点。
这意味着：原本需要为 8 个点开工的子工厂，现在只需要为 4 个点开工。

3. 对称性（Symmetry Property） —— 蝴蝶碰撞的由法

这是你之前看到的“一加一减”逻辑的物理来源。

数学表示：$\omega_n^{k + n/2} = -\omega_n^k$
演示：
在圆周上，走半圈（$n/2$ 步）到达的点，恰好是起点的镜像点（相反数）。
- 例如 $N=8$ 时，$\omega_8^4 = -\omega_8^0 = -1$。
- 例如 $N=8$ 时，$\omega_8^5 = -\omega_8^1$。
工业意义：
这保证了我们在计算 $A(x) = \text{偶项} + x \cdot \text{奇项}$ 时，对于对称的点只需要计算一次乘法。加法给正点，减法给负点。

4. 消去性（Summation Property / Orthogonality） —— 逆变换的保障

这是最神奇的性质，它保证了我们能把“乱七八糟”的频域数字变回 56088。

数学表示：
\[\sum_{j=0}^{n-1} (\omega_n^k)^j = \begin{cases} n, & \text{if } k \equiv 0 \pmod n \\ 0, & \text{otherwise} \end{cases}\]
直观理解：
如果你把所有的单位根加起来，它们会因为受力平衡而抵消，结果是 $0$。只有当所有探测点“步调一致”（全部是 1）时，结果才会叠加。
工业应用：
这个性质被用于逆变换（INTT）。它就像一把特殊的筛子，能够从混合的频率信号中，精确地“筛”出原始的系数 $a_k$。这也是为什么逆变换只需要把 $\omega$ 换成 $\omega^{-1}$ 并在最后除以 $N$ 就能还原数据的理论根据。

  
const nttMultiply = (num1, num2) => {
    const P = 998244353n;
    const G = 3n; // 模 P 的原根

    // 快速幂运算：计算 (base^exp) % P
    const power = (a, b) => {
        let res = 1n;
        a %= P;
        while (b > 0n) {
            if (b % 2n === 1n) res = (res * a) % P;
            a = (a * a) % P;
            b /= 2n;
        }
        return res;
    };

    // 计算乘法逆元：利用费马小定理 a^(P-2) % P
    const modInverse = (n) => power(n, P - 2n);

    const a = num1.toString().split('').map(x => BigInt(x)).reverse();
    const b = num2.toString().split('').map(x => BigInt(x)).reverse();

    let n = 1;
    while (n < a.length + b.length) n <<= 1;

    // 这个地方已经把长度补到二次幂了 并且是一定能容纳乘积的
    // a位数乘b位数最多就是(a+b)位
    const fa = new BigInt64Array(n);
    const fb = new BigInt64Array(n);
    a.forEach((v, i) => fa[i] = v);
    b.forEach((v, i) => fb[i] = v);

    // 核心 NTT 变换：逻辑与 FFT 完全相同
    const ntt = (a, invert) => {
        const n = a.length;
        // 位反转排序 (Bit-reversal permutation)
        // 这个排序的效果是没有直观的解释的 唯一目的就是方便下面的蝴蝶操作
        for (let i = 1, j = 0; i < n; i++) {
            let bit = n >> 1;
            for (; j & bit; bit >>= 1) j ^= bit;
            j ^= bit;
            if (i < j) [a[i], a[j]] = [a[j], a[i]];
        }

        // 蝴蝶操作
        // len is 2 4 8 16 ...
        for (let len = 2; len <= n; len <<= 1) {
            // 这里用原根代替了复数域的 e^(-i*2*pi/n)
            // wlen 其实是一个比较固定的数字序列 因为GP都是固定的 然后len是2的次幂
            let wlen = power(G, (P - 1n) / BigInt(len));
            if (invert) wlen = modInverse(wlen);

            // i is 0 2 4 8 ... length of one of the input number
            for (let i = 0; i < n; i += len) {
                let w = 1n;
                // j is 0 1 2 3 ... 2^(some number)-1
                for (let j = 0; j < len / 2; j++) {
                    // one of the coeffs of our polynominal
                    const u = a[i + j];
                    const v = (a[i + j + len / 2] * w) % P;
                    a[i + j] = (u + v) % P;
                    a[i + j + len / 2] = (u - v + P) % P;
                    w = (w * wlen) % P;
                }
            }
        }

        if (invert) {
            const nInv = modInverse(BigInt(n));
            for (let i = 0; i < n; i++) a[i] = (a[i] * nInv) % P;
        }
    };

    ntt(fa, false);
    ntt(fb, false);

    // 频域点乘：逐对相乘 (Hadamard Product)
    const fc = new BigInt64Array(n);
    for (let i = 0; i < n; i++) fc[i] = (fa[i] * fb[i]) % P;

    ntt(fc, true);

    // 进位处理
    const result = [];
    let carry = 0n;
    for (let i = 0; i < n; i++) {
        let val = fc[i] + carry;
        result.push(Number(val % 10n));
        carry = val / 10n;
    }
    while (carry) {
        result.push(Number(carry % 10n));
        carry /= 10n;
    }

    while (result.length > 1 && result[result.length - 1] === 0) result.pop();
    return result.reverse().join('');
};

console.log(nttMultiply("506", "208")); // 输出: 105248

本文由作者按照 CC BY 4.0 进行授权