FFT解大整数乘法

2018-12-13

1. 代码文件说明

提交的代码包括两个文件夹，fft 和 judge。fft 文件夹包括的是可独立运行的 fft 求大数乘法的程序（做了输入输出重定向），以及实际提交于 Leetcode 的代码。judge 文件夹包含了三个数据文件（test.data 为十组本地测试数据，ans.data 为正确答案，fftans.data 为 fft.exe 的运行输出结果），两个 python 文件（data.py 为生成测试数据和答案的脚本，judge.py 为判断结果是否正确的脚本）。

代码放在文末。

2. 求解思想

系数表示法

在计算乘法的时候，我们可以把一个数字分解为一个多项式 $A(x) = \sum^{n-1}_{i=0}\alpha_ix^i$ 。对于一个固定的 $x$，我们就可以把两个数字分解为 $n$, $m$ 维的两个向量（两个 $x$ 进制数），则乘法结果就是对应的 $n$ 维和 $m$ 维的向量的卷积结果（多项式系数的卷积）。根据 卷积定理 ：

向量卷积的离散傅里叶变换是向量离散傅里叶变换的乘积。

我们就可以把两个大数分解为向量 $a$，$b$，然后对两个向量分别做离散傅里叶变换，之后逐位相乘得到一个新向量 $c$。对 $c$ 做逆离散傅里叶变换，再进行进位，就可以得到结果。

快速傅里叶变换可以做到在 $O(NlogN)$ 的复杂度完成原本 $O(N^2)$ 复杂度的离散傅里叶变换，让整个算法的复杂度也降低到 $O(NlogN)$ 。
点值表示法

对 卷积定理 进一步进行展开解释。在两个整数分解为 $A(x)$ 的时候，我们可以分别对于 $n$ 个不同的 $x$ 值记录 $n$ 个点值对 $(x_i,A_a(x_i))$ 和 $(x_i,A_b(x_i))$ 。那么我们也可以根据这两组点值对直接得到乘法结果对应的一组点值对 $(x_i,A_a(x_i)\times A_b(x_i))$ 。

多项式插值的唯一性定理 证明了这 $n$ 个点值对即可用来表示（还原）对应的多项式。由于乘法结果会有 $2n$ 位，所以我们一开始分解时实际分别记录两个数对应的 $2n$ 个点值对，最终可以使用插值函数的方法恢复出 $2n$ 位的乘法结果。

这和之前系数表示法那一段之间的联系就在于，离散傅里叶变换 DFT 就是一个计算系数表达式与点值表达式之间互换的算法（采样）。只不过实际计算过程中，DFT 选取的 $x$ 值是复数值（n次单位根）。

直接使用 DFT 计算的话，转换步骤的复杂度显然需要 $O(n^2)$ 。而 FFT 能够快速完成系数表达式和点值表达式之间的换算，使复杂度降低到 $O(NlogN)$ 。
FFT 计算过程

结合 消去引理，折半引理，求和引理，我们可以分治的求解 DFT：

一个界长为 $N$ 的离散傅里叶变换可以重新写成两个界长各为 $N/2$ 的离散傅里叶变换之和。其中一个变换由原来 $N$ 个点中的偶数点构成，另一个变换由奇数点构成。这个过程递归进行下去，直到将全部数据细分为界长为 1 的变换。在边界上，界长为 1 的变换等于自身。

为了方便起见，在反转置换之前先补充前导零，使向量长度为 2 的幂形式。然后对向量进行裂项，即反转变换。将一个 2 的 n 次幂长度的向量进行裂项操作，每个元素的位置就会是下标的二进制反转之后再转换成十进制的位置，可以采用 Rader 算法（二进制平摊反转置换算法）实现，这样可以将递归转化成迭代执行。

最终在回溯的时候，不断套用公式执行蝶形运算，就可以得到结果了。

由于递归执行了求解过程，递归树为一棵完全二叉树，所以易知复杂度为 $O(NlogN)$ 。
逆 FFT

只需对原来的 FFT 算法代码进行小小的修改：a 和 y 互换，$ω_n^{-1}$ 代替 $ω_n$ ，最后所有值除以 $n$ 即可。
全部程序执行流程：

读入两个大数字符串，分别转化为数组一位一位存储，并补充前导零。分别进行 FFT 计算，按位相乘，然后做 IFFT 计算。最后将结果进行进位操作，就能得到最终的乘法结果。

3. 代码执行结果

本地数据

共生成了 10 组 100 $\times$ 100 位的随机数。

示例：（第一组大数）

7739385993211797423647071118580282469713569881037743170530795280641276969768173826242862186300508114

9672516198036485560430536046045403561144663114397844686576323489397779756322778671971277864423561283

得到结果：

74859336382197804460271536901670667372336115116031816158713680432079139311256687277234225158377643751548229859235061475750266007841281857699781669159860971180777171200406847628409854302216736317750262

使用 Python 脚本判断十组数据运行正确性：

通过测试。
Leetcode (https://leetcode.com/problems/multiply-strings/)

通过。

P.S. Leetcode 上这道题最快的 4 ms 过题代码使用的是最朴素的 $O(n^2)$ 逐位相乘然后进位的算法。我的代码慢了 50 多倍，这说明 FFT 的常数大概是非常大了（至少我这个代码是）。

代码

fft/fft.cpp

  1#include <cstdio>
  2#include <iostream>
  3#include <cmath>
  4#include <cstring>
  5#include <algorithm>
  6
  7using namespace std;
  8
  9const double PI = acos(-1.0);
 10
 11// 复数结构体
 12struct Complex{
 13    double x, y;     // 实部和虚部 x + yi
 14    Complex(double _x = 0.0, double _y = 0.0){
 15        x = _x;
 16        y = _y;
 17    }
 18    Complex operator - (const Complex &b) const{
 19        return Complex(x - b.x, y - b.y);
 20    }
 21    Complex operator + (const Complex &b) const{
 22        return Complex(x + b.x, y + b.y);
 23    }
 24    Complex operator * (const Complex &b) const{
 25        return Complex(x*b.x - y*b.y, x*b.y + y*b.x);
 26    }
 27};
 28
 29// FFT 和 IFFT 前的反转变换。
 30// 位置 i 和 （i 二进制反转后位置）互换
 31// len 必须是 2 的幂
 32void change(Complex y[], int len){
 33    int i, j, k;
 34    for(i = 1, j = len / 2; i < len - 1; i++){
 35        if(i < j) swap(y[i], y[j]);
 36        // 交换互为小标反转的元素，i < j 保证交换一次
 37        // i 做正常的 +1，j 左反转类型的 +1,始终保持 i 和 j 是反转的
 38        k = len / 2;
 39        while(j >= k){
 40            j -= k;
 41            k /= 2;
 42        }
 43        if(j < k) j += k;
 44    }
 45}
 46
 47// FFT
 48// len 必须为 2 的幂，
 49// on == 1 时是 DFT，on == -1 时是 IDFT
 50void fft(Complex y[], int len, int on){
 51    change(y, len);
 52    for(int h = 2; h <= len; h *= 2){
 53        Complex wn(cos(-on * 2*PI / h), sin(-on * 2*PI / h));
 54        for(int j = 0; j < len; j += h){
 55            Complex w(1, 0);
 56            for(int k = j; k < j + h/2; k++){
 57                Complex u = y[k];
 58                Complex t = w * y[k + h/2];
 59                y[k] = u + t;
 60                y[k + h/2] = u - t;
 61                w = w * wn;
 62            }
 63        }
 64    }
 65    if(on == -1)
 66        for(int i = 0; i < len; i++)
 67            y[i].x /= len;
 68}
 69
 70const int MAXN = 200010;
 71Complex x1[MAXN], x2[MAXN];
 72char str1[MAXN / 2], str2[MAXN / 2];
 73int sum[MAXN];
 74
 75int main(){
 76    freopen("test.data", "r", stdin);
 77    freopen("fftans.data", "w", stdout);
 78
 79    // 输入两个大数字符串
 80    while(scanf("%s%s", str1, str2)!=EOF){
 81        int len1 = strlen(str1);
 82        int len2 = strlen(str2);
 83
 84        int len = 1;
 85        while(len < len1 * 2 || len < len2 * 2) len *= 2;
 86
 87        // 把字符串转化为复数数组并补前导 0 至 2 的幂
 88        for(int i = 0; i < len1; i++)
 89            x1[i] = Complex(str1[len1 - 1 - i] - '0', 0);
 90        for(int i = len1; i < len; i++)
 91            x1[i] = Complex(0, 0);
 92
 93        for(int i = 0; i < len2; i++)
 94            x2[i] = Complex(str2[len2 - 1 - i] - '0', 0);
 95        for(int i = len2; i < len; i++)
 96            x2[i] = Complex(0, 0);
 97        
 98        // 求DFT
 99        fft(x1, len, 1);
100        fft(x2, len, 1);
101
102        // 相乘
103        for(int i = 0; i < len; i++)
104            x1[i] = x1[i] * x2[i];
105        
106        // IDFT
107        fft(x1, len, -1);
108
109        // 转化为整数值
110        for(int i = 0; i < len; i++)
111            sum[i] = (int)(x1[i].x + 0.5);
112        for(int i = 0; i < len; i++){
113            sum[i+1] += sum[i] / 10;
114            sum[i] %= 10;
115        }
116
117        // 输出结果
118        len = len1 + len2 - 1;
119        while(sum[len] <= 0 && len > 0) len--;
120        for(int i = len; i >= 0; i--)
121            printf("%c", sum[i] + '0');
122        printf("\n");
123    }
124    return 0;
125}

fft/leetcode_upload.cpp

  1class Solution {
  2public:
  3    double PI = acos(-1.0);
  4    // 复数结构体
  5    struct Complex{
  6        double x, y;     // 实部和虚部 x + yi
  7        Complex(double _x = 0.0, double _y = 0.0){
  8            x = _x;
  9            y = _y;
 10        }
 11        Complex operator - (const Complex &b) const{
 12            return Complex(x - b.x, y - b.y);
 13        }
 14        Complex operator + (const Complex &b) const{
 15            return Complex(x + b.x, y + b.y);
 16        }
 17        Complex operator * (const Complex &b) const{
 18            return Complex(x*b.x - y*b.y, x*b.y + y*b.x);
 19        }
 20    };
 21    // FFT 和 IFFT 前的反转变换。
 22    // 位置 i 和 （i 二进制反转后位置）互换
 23    // len 必须是 2 的幂
 24    void change(Complex y[], int len){
 25        int i, j, k;
 26        for(i = 1, j = len / 2; i < len - 1; i++){
 27            if(i < j) swap(y[i], y[j]);
 28            // 交换互为小标反转的元素，i < j 保证交换一次
 29            // i 做正常的 +1，j 左反转类型的 +1,始终保持 i 和 j 是反转的
 30            k = len / 2;
 31            while(j >= k){
 32                j -= k;
 33                k /= 2;
 34            }
 35            if(j < k) j += k;
 36        }
 37    }
 38
 39    // FFT
 40    // len 必须为 2 的幂，
 41    // on == 1 时是 DFT，on == -1 时是 IDFT
 42    void fft(Complex y[], int len, int on){
 43        change(y, len);
 44        for(int h = 2; h <= len; h *= 2){
 45            Complex wn(cos(-on * 2*PI / h), sin(-on * 2*PI / h));
 46            for(int j = 0; j < len; j += h){
 47                Complex w(1, 0);
 48                for(int k = j; k < j + h/2; k++){
 49                    Complex u = y[k];
 50                    Complex t = w * y[k + h/2];
 51                    y[k] = u + t;
 52                    y[k + h/2] = u - t;
 53                    w = w * wn;
 54                }
 55            }
 56        }
 57        if(on == -1)
 58            for(int i = 0; i < len; i++)
 59                y[i].x /= len;
 60    }
 61
 62    Complex x1[200010], x2[200010];
 63    char str1[100005], str2[100005];
 64    int sum[200010];
 65
 66    string multiply(string num1, string num2) {
 67        // 输入两个大数字符串
 68        strcpy(str1, num1.c_str());
 69        strcpy(str2, num2.c_str());
 70
 71        int len1 = strlen(str1);
 72        int len2 = strlen(str2);
 73
 74        int len = 1;
 75        while(len < len1 * 2 || len < len2 * 2) len *= 2;
 76
 77        // 把字符串转化为复数数组
 78        for(int i = 0; i < len1; i++)
 79            x1[i] = Complex(str1[len1 - 1 - i] - '0', 0);
 80        for(int i = len1; i < len; i++)
 81            x1[i] = Complex(0, 0);
 82        for(int i = 0; i < len2; i++)
 83            x2[i] = Complex(str2[len2 - 1 - i] - '0', 0);
 84        for(int i = len2; i < len; i++)
 85            x2[i] = Complex(0, 0);
 86        
 87        // 求DFT
 88        fft(x1, len, 1);
 89        fft(x2, len, 1);
 90
 91        // 相乘
 92        for(int i = 0; i < len; i++)
 93            x1[i] = x1[i] * x2[i];
 94        
 95        // IDFT
 96        fft(x1, len, -1);
 97
 98        // 转化为整数值
 99        for(int i = 0; i < len; i++)
100            sum[i] = (int)(x1[i].x + 0.5);
101        for(int i = 0; i < len; i++){
102            sum[i+1] += sum[i] / 10;
103            sum[i] %= 10;
104        }
105
106        // 输出结果
107        len = len1 + len2 - 1;
108        while(sum[len] <= 0 && len > 0) len--;
109        string ans;
110        for(int i = len; i >= 0; i--)
111            ans.append(1, char(sum[i] + '0'));
112        return ans;
113    }
114};

judge/data.py

 1from random import randint
 2
 3test = open("test.data", 'w')
 4ans = open("ans.data", 'w')
 5
 6last = str()
 7temp = str()
 8
 9for i in range(20):
10    last = temp
11    temp = str(randint(1, 9))
12    for j in range(99):
13        temp += str(randint(0, 9))
14    test.write(temp)
15    if i % 2 == 1:
16        test.write("\n")
17        ans.write(str(int(temp) * int(last)) + "\n")
18    else:
19        test.write(" ")
20
21test.close()

judge/judge.py

 1fftans = open("fftans.data", "r").readlines()
 2ans = open("ans.data", "r").readlines()
 3
 4if len(fftans) != len(ans):
 5    print("Fail")
 6
 7for i in range(len(fftans)):
 8    if fftans[i] != ans[i]:
 9        print("Fail ! No. {}".format(i + 1))
10        break
11print("Success !")

#courses
本文总字数 2990
本文阅读量
本站访客量

↪ reply by email