【fft】使用快速傅里叶变换计算大整数乘法

使用快速傅里叶变换计算大整数乘法

这是一道网络教室的题目，本来是要用分块乘法解决的，但是快速傅里叶也是一种解法。

1. 从多项式入手

实数可以表示成为多项式的形式，举一个例子
$$143=3\times 10^0+4\times 10^1+1\times 10^2$$

这说明 $143$ 是多项式
$$P(x)=3x^0+4x+1x^2$$

是当 $x=10$ 时的特殊情况。更一般地，所有的实数乘法都可以看作多项式相乘的形式，解决了多项式相乘的问题就解决了数乘的问题。这看似将问题复杂化了，但是不要着急，看下去就会有结果。

有两个多项式
$$\begin{gathered} A(x)=\sum _{i=0}^{s-1}{a}_i{x}^i \\ B(x)=\sum _{i=0}^{t-1}{b}_i{x}^i \end{gathered}$$

按照分配律计算 $A\times B$ 则
$$\begin{gathered} C(x)=\sum _{i=0}^{t+s-2}{c}_i{x}^i \\ {c}_i=\sum _{j=0}^i\sum _{k=0}^{i-j}{a}_j{b}_k \end{gathered}$$

大致的时间复杂度为 $\Theta (st)$。

2. 换一种方式看待多项式

书接上文，多项式 $A$ 的计算结果是由输入的变量 $x$ 和系数 $a_0,\cdots ,a_{s-1}$ 决定的，所以在表达多项式 $A$ 时可以略去输入 $x$，从而看作 $s$ 个系数。对于不同的 $s$ 个输入 $x_0,\cdots ,x_{s-1}$ 我们可以写出
$$y_i=A(x_i)=a_0+a_1{x}_i+\cdots +a_{s-1}{x}_i^{s-1},i=0,\cdots ,s-1$$

这是可以写成矩阵形式的
$$\left[\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^{s-1}\\ 1& x_1& x_1^2& \cdots & x_1^{s-1}\\ & & ⋮& & \\ 1& x_{s-1}& x_{s-1}^2& \cdots & x_{s-1}^{s-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{s-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{s-1}\end{array}\right]$$

由于范德蒙行列式
$$\left|\begin{array}{ccccc}1& x_0& x_0^2& \cdots & x_0^{s-1}\\ 1& x_1& x_1^2& \cdots & x_1^{s-1}\\ & & ⋮& & \\ 1& x_{s-1}& x_{s-1}^2& \cdots & x_{s-1}^{s-1}\end{array}\right|\ne 0$$

故对于多项式 $A$ 即使在不知道系数的取值但是知道 $x_i,y_i$ 的取值的情况下也可以通过矩阵逆运算得出系数的值，$a=V^{-1}y$。所以多项式 $A$ 也可以表达为自变量取值-多项式值对的形式，称为点值表达。
A → { ( x 0 , y 0 ) , ⋯   , ( x s − 1 , y s − 1 ) } B → { ( x 0 , y 0 ′ ) , ⋯   , ( x t − 1 , y t − 1 ′ ) } A\to\{(x_0,y_0),\cdots,(x_{s-1},y_{s-1})\}\\ B\to\{(x_0,y_0'),\cdots,(x_{t-1},y_{t-1}')\} A→{(x0​,y0​),⋯,(xs−1​,ys−1​)}B→{(x0​,y0′​),⋯,(xt−1​,yt−1′​)}

这样表达带来了一个优势，那就是乘法的速度变快了。我们假设 $s=t=m$（即使 $s\ne t$也可以通过添加零系数的方式添到一样），那么
C → { ( x 0 , y 0 y 0 ′ ) , ⋯   , ( x n − 1 , y m − 1 y m − 1 ′ ) } C\to\{(x_0,y_0y_0'),\cdots,(x_{n-1},y_{m-1}y_{m-1}')\} C→{(x0​,y0​y0′​),⋯,(xn−1​,ym−1​ym−1′​)}

的计算时间复杂度在知道 $A,B$ 的点值表达的情况下为 $\Theta (n)$。但是我们需要注意对于多项式 $C$ 拥有最多 $2m$ 项系数，所以如果需要通过矩阵求逆 $c=V^{-1}y$ 的方式求出多项式 $C$ 的常规表达法，则需要 $2m$ 个不同的 $x$ 取值。

那么现在任务就明确了，对于两个输入的大整数（A, B）

1. 我们按照位数写成多项式的形式，然后求出对应且足够的点值表达
2. 然后计算出相乘后的点值表达（C），时间复杂度为 $\Theta (m)$
3. 最后通过矩阵求逆的运算计算出点值表达对应的系数
4. 对系数进行整理从而得到整数乘法的结果。

3. 单位复根

为了加速运算，我们引入单位复根这一概念。它的使用将在后面看到。

n次单位复根：使得 $w^n=1,n\in N$ 的复数 $w$.

主n次单位复根：$w_n=\exp (2\pi i/n)$.

这里给出几个简单的性质以供使用和参考：

消去定律：$w_{dn}^{dk}=w_n^k.$

折半定律：$w_n^{n/2}=w_2=-1,n\text{ is even.}$

求和定律：对于任意整数 $n⩾1$ 和非负整数 $k,\gcd (k,n)=1$，有
$$\sum _{j=0}^{n-1}(w_n^k{)}^j=0.$$
求和定律的证明可以使用等比数列求和得到。

4. 求点值表达

对于多项式 $A,B$，我们首先通过添加零系数的方式将它们的项数补充到一致，然后再通过添加零系数的方式将两个多项式的项数都添加到一个 $2$ 的次幂的数上，将这个数翻倍，设为 $n$。这里翻倍的原因是为了使得方程有解，原因在前文中已经说过。

令 $x_i=w_n^i$，我们要求出
$$\begin{gathered} A(w_n^0),A(w_n^1),\cdots ,A(w_n^{n-1}) \\ B(w_n^0),B(w_n^1),\cdots ,B(w_n^{n-1}) \end{gathered}$$

以下以 $A(x)$ 为例，由于扩充后 $A$ 的项数为 $2$ 的次幂，我们设
$$\begin{gathered} A^{[0]}=a_0+a_{2}x+\cdots +a_{n-2}{x}^{n/2-1} \\ {A}^{[1]}=a_1+a_{3}x+\cdots +a_{n-1}{x}^{n/2-1} \end{gathered}$$

则有 $A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)$，给出其计算复杂度 $T(n)=2T(n/2)+\Theta (n)=\Theta (n\lg n)$。

特别说明一下，设 $y_k,y_k^{[0]},y_k^{[1]},k=0,\cdots ,n/2-1$ 为多项式 $A(w_n^k),A^{[0]}(w_n^{2k}),A^{[1]}(w_n^{2k})$ 的值，那么
$$\begin{array}{rl}{y}_k& =y_k^{[0]}+w_n^k{y}_k^{[0]}\\ & =A^{[0]}(w_n^{2k})+w_n^k{A}^{[1]}(w_n^{2k})\\ & =A(w_n^k)\\ y_{k+n/2}& =y_k^{[0]}-w_n^k{y}_k^{[0]}\\ & =y_k^{[0]}+w_n^{k+n/2}{y}_k^{[0]}\\ & =A^{[0]}(w_n^{2k+n})+w_n^{k+n/2}{A}^{[1]}(w_n^{2k+n})\\ & =A^{[0]}(w_n^{2k})-w_n^k{A}^{[1]}(w_n^{2k})\\ & =A(w_n^{k+n/2})\end{array}$$

这指导我们我们只需遍历 $k=0,1,\cdots ,n/2-1$ 就可以计算出 $A(w_n^0),\cdots ,A(w_n^{n-1})$。

5. 求系数

通过上述方式我们可以求出多项式 $C=A\times B$ 的点表达式，且时间复杂度为 $\Theta (n\lg n)$，接下来需要通过点表达式求出各项系数，这关键是要求出逆矩阵乘法的结果，即
$$\left[\begin{array}{c}{c}_0\\ ⋮\\ c_{n-1}\end{array}\right]=V(w_n^0,w_n^1,\cdots ,w_n^{n-1}{)}^{-1}\left[\begin{array}{c}{y}_0\\ ⋮\\ y_{n-1}\end{array}\right]$$

由于单位复根的性质，可以证明 $V(w_n^0,w_n^1,\cdots ,w_n^{n-1}{)}^{-1}=[w_n^{-jk}/n{]}_{n\times n}$，即这个范德蒙矩阵的逆矩阵的第 $j$ 行第 $k$ 列的元素为 $w_n^{-jk}/n$。

相当于
$$\left[\begin{array}{c}{c}_0\\ ⋮\\ c_{n-1}\end{array}\right]=\frac{1}{n}V(w_n^0,w_n^{-1},\cdots ,w_n^{-(n-1)})\left[\begin{array}{c}{y}_0\\ ⋮\\ y_{n-1}\end{array}\right]$$

利用之前的方法同样可以在 $\Theta (n\lg n)$ 的时间内完成运算。

6. 例子

这里举例说明以上的运算方式
$$\begin{gathered} a=123,b=89 \\ a\times b \end{gathered}$$

构建出两个多项式
$$A(x)=3+2x+x^2,B=9+8x$$

根据位数，最终的结果的位数有可能达到 $3+2=7$ 位，向上填充到 $2$ 的次幂数 $8$，所以我们需要8次单位复根。
$$\exp (i\pi /4\times 0),\exp (i\pi /4\times 1),\cdots ,\exp (i\pi /4\times 7)$$

分别带入 $A,B$ 得到 $C$

w	A	B	C
$w_8^0$	$6.000+0.000i$	$17.000+0.000i$	$102.000+0.000i$
$w_8^1$	$4.414+2.414i$	$14.657+5.657i$	$51.042+60.355i$
$w_8^2$	$2.000+2.000i$	$9.000+8.000i$	$2.000+34.000i$
$w_8^3$	$1.586+0.414i$	$3.343+5.657i$	$2.958+10.355i$
$w_8^4$	$2.000+0.000i$	$1.000+0.000i$	$2.000+0.000i$
$w_8^5$	$1.586-0.414i$	$3.343-5.657i$	$2.958-10.355i$
$w_8^6$	$2.000-2.000i$	$9.000-8.000i$	$2.000-34.000i$
$w_8^7$	$4.414-2.414i$	$14.657-5.657i$	$51.042-60.355i$

所以对于多项式 $C$ 的点值表达有
$$\left[\begin{array}{cccccccc}1.000+0.000i& 1.000+0.000i& 1.000+0.000i& 1.000+0.000i& 1.000+0.000i& 1.000+0.000i& 1.000+0.000i& 1.000+0.000i\\ 1.000+0.000i& 0.707+0.707i& 0.000+1.000i& -0.707+0.707i& -1.000+0.000i& -0.707-0.707i& -0.000-1.000i& 0.707-0.707i\\ 1.000+0.000i& 0.000+1.000i& -1.000+0.000i& -0.000-1.000i& 1.000+0.000i& 0.000+1.000i& -1.000+0.000i& -0.000-1.000i\\ 1.000+0.000i& -0.707+0.707i& -0.000-1.000i& 0.707+0.707i& -1.000+0.000i& 0.707-0.707i& 0.000+1.000i& -0.707-0.707i\\ 1.000+0.000i& -1.000+0.000i& 1.000+0.000i& -1.000+0.000i& 1.000+0.000i& -1.000+0.000i& 1.000+0.000i& -1.000+0.000i\\ 1.000+0.000i& -0.707-0.707i& 0.000+1.000i& 0.707-0.707i& -1.000+0.000i& 0.707+0.707i& -0.000-1.000i& -0.707+0.707i\\ 1.000+0.000i& -0.000-1.000i& -1.000+0.000i& 0.000+1.000i& 1.000+0.000i& -0.000-1.000i& -1.000+0.000i& 0.000+1.000i\\ 1.000+0.000i& 0.707-0.707i& -0.000-1.000i& -0.707-0.707i& -1.000+0.000i& -0.707+0.707i& 0.000+1.000i& 0.707+0.707i\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\\ c_7\end{array}\right]=\left[\begin{array}{c}102.000+0.000i\\ 51.042+60.355i\\ 2.000+34.000i\\ 2.958+10.355i\\ 2.000+0.000i\\ 2.958-10.355i\\ 2.000-34.000i\\ 51.042-60.355i\end{array}\right]$$

求解得到
$$\left[\begin{array}{c}{c}_0\\ c_1\\ c_2\\ c_3\\ c_4\\ c_5\\ c_6\\ c_7\end{array}\right]=\left[\begin{array}{c}27\\ 42\\ 25\\ 8\\ 0\\ 0\\ 0\\ 0\end{array}\right]$$

即 $C(x)=27+42x+25x^2+8x^3=A(x)\times B(x)\Rightarrow C(10)=A(10)\times B(10)\Rightarrow 27+42\times 10+25\times 10^2+8\times 10^3=10947=123\times 89$.

7. 代码

#include <iostream>
#include <complex>
#include <cmath>
#include <vector>

using namespace std;
typedef complex<double> dcomp;

class mypoly {
public:
    size_t len;
    int *memory;

    explicit mypoly(size_t l) : memory(new int[l]), len(l) {}

    mypoly(mypoly &&x) noexcept {
        len = x.len;
        memory = x.memory;
        x.memory = nullptr;
    }

    ~mypoly() {
        delete[] memory;
    }

    void read(string str) const {
        for (size_t i = 0; i < str.length(); i++)
            memory[i] = str[str.length() - i - 1] - '0';
        for (size_t i = str.length(); i < len; i++)
            memory[i] = 0;
    }
};

// 正向傅里叶变换
vector<dcomp> fft(const mypoly &a, size_t off, size_t sep) {
    if (sep == a.len) {
        dcomp a0 = a.memory[off];
        vector<dcomp> r;
        r.push_back(a0);
        return r;
    } else {
        size_t n = a.len / sep;
        dcomp wn(cos(2 * M_PI / n), sin(2 * M_PI / n)), w(1., 0.);
        vector<dcomp> y0 = fft(a, off, sep << 1)
                , y1 = fft(a, off + sep, sep << 1)
                        , y(a.len / sep);
        for (size_t k = 0; k < y.size() >> 1; k++) {
            y[k] = y0[k] + w * y1[k];
            y[k + (n >> 1)] = y0[k] - w * y1[k];
            w = w * wn;
        }
        return y;
    }
}

// 反向傅里叶变换
vector<dcomp> ifft(const vector<dcomp> &y, size_t off, size_t sep) {
    if (sep == y.size()) {
        dcomp y0 = y[off];
        vector<dcomp> r;
        r.push_back(y0);
        return r;
    } else {
        size_t n = y.size() / sep;
        dcomp wn(cos(2 * M_PI / n), -sin(2 * M_PI / n)), w(1., 0.);
        vector<dcomp> a0 = ifft(y, off, sep << 1)
                , a1 = ifft(y, off + sep, sep << 1)
                        , a(y.size() / sep);
        for (size_t k = 0; k < a.size() >> 1; k++) {
            a[k] = a0[k] + w * a1[k];
            a[k + (n >> 1)] = a0[k] - w * a1[k];
            w = w * wn;
        }
        return a;
    }
}

int main() {
    string str1, str2;
    cin >> str1 >> str2;
    size_t m = str1.size() | str2.size(), l = 1;
    while (l <= m) l <<= 1;
    l <<= 1;
    mypoly p1(l), p2(l);
    p1.read(str1), p2.read(str2);

    vector<dcomp> r1 = fft(p1, 0, 1), r2 = fft(p2, 0, 1), r(l);
    for (size_t ind = 0; ind < l; ind++)
        r[ind] = r1[ind] * r2[ind];

    vector<dcomp> a = ifft(r, 0, 1);
    vector<int> out(a.size());
    for (size_t k = 0; k < out.size(); k++)
    	// 因为整数截断的原因所以加0.5再向下取整防止诸如1.9998被截断成1而不是2
        out[k] = (int)(std::abs(a[k]) / l + .5);
    for (size_t k = 0; k < out.size() - 1; k++) {
        out[k + 1] += out[k] / 10;
        out[k] %= 10;
    }
    size_t end = out.size();
    while(out[--end] == 0);
    while(end != 0)
        cout << out[end--];
    cout << out[0] << endl;

    return 0;
}

匆匆写就，如有错误恳请指出。