基于FFT的大整数乘法

最新推荐文章于 2023-04-05 16:06:55 发布

gokingd

最新推荐文章于 2023-04-05 16:06:55 发布

阅读量2.2k

点赞数 6

分类专栏：算法文章标签：算法算法导论 fft

本文链接：https://blog.csdn.net/qq_43290092/article/details/121043417

版权

算法专栏收录该内容

19 篇文章 0 订阅

订阅专栏

本文介绍了快速傅里叶变换(FFT)的基本原理及其在计算多项式值和求解多项式系数中的高效性。通过矩阵形式展示了如何从已知点的值反推出多项式系数，并详细解释了如何通过FFT加速多项式乘法，特别是在计算大整数乘法时的应用。同时，提供了相关LeetCode题目的解题思路和代码实现，强调了FFT在处理这类问题时的分治算法本质。

摘要由CSDN通过智能技术生成

多项式求值
对于多项式 $f(x)=a_{0}+a_{1}*x^{1}+a_{2}*x^{2}+...+a_{n-1}*x^{n-1}$ 通常我们要花费 $O(n^2)$ 的时间来求其 $n$ 个点的处的取值。1965年Cooley and Tukey提出的FFT能够在 $O (n l o g n)$ 的时间求得 $n$ 个点处的值。其本质是利用复数的周期性和分治算法消除计算的冗余。故FFT加速求值时，只能计算特殊点处的取值，如 $1,w,w^2,...,w^{n-1}$ ,其中 $w=e^{2\pi*i/n}$
具体过程这里不做推导。直接给出FFT计算多项式值的伪代码。

--引用卜算的PPT

–引用卜算的PPT

求多项式
那么如果已知 $n$ 个点的值，如何求 $f (x)$ 呢？首先给出求多项式的值的矩阵形式。（~~本来想用Latex写，偷懒一下，直接截图算了~~ ）
在这里插入图片描述

通过矩阵形式，显然可以看出，我们在已知 $n$ 个点的值时，通过求矩阵的逆左乘多项式的值，就可以得到多项式的系数。通常情况下，采用高斯消元求矩阵的逆，需要 $O(n^3)$ 的时间。但是如果是特定点处的指，就可以直接通过复数的性质得到逆矩阵，同时利用复数的周期性和分治算法，类似于FFT，就可以在 $O (n l o g n)$ 的时间求得多项式的系数。

在这里插入图片描述

所以只需要在FFT的代码上稍加修改即可。

多项式乘法
对于多项式 $A (x)$ 和 $B (x)$ 求 $C (x)$
$A(x)=a_{0}+a_{1}*x^{1}+a_{2}*x^{2}+...+a_{n-1}*x^{n-1} \\ B(x)=b_{0}+b_{1}*x^{1}+b_{2}*x^{2}+...+b_{n-1}*x^{n-1} \\ C(x)=A(x)*B(x)=c_{0}+c_{1}*c^{1}+c_{2}*x^{2}+...+c_{2n-2}*x^{2n-2}$
$c_{k}=\Sigma_{i=0}^{k}a_{i}b_{k-i}$ 的每一项其实是计算卷积，求出C需要的时间为 $O(n^2)$ 。这时候就要放上一张经典图了。
在这里插入图片描述
从我们上面的逻辑来看，就是先通过FFT计算处特定点处的值，得到多项式C在特定点处的值，再通过IFFT求得C的系数。其实本质上，A*B求C的过程是求卷积的过程，而我们可以通过FFT来加速求卷积的过程。而图上的过程，其实是通过FFT，转换到频域下进行计算，FFT其实在信号领域，就是求其频域的值，只不过在目前的算法里面，暂时不需要去想时域和频域。
给一道用FFT求多项式乘积的题。洛谷P3803，FFT模板题。传送门

#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
const int MAX = 4e6 + 10;
struct Complex{
    double x, y;
    Complex(){x=y=0;}
    Complex(double x,double y){
        this->x = x;
        this->y = y;
    }
}a[MAX], b[MAX];
Complex operator+(Complex a, Complex b) { return Complex(a.x + b.x, a.y + b.y); }
Complex operator-(Complex a, Complex b) { return Complex(a.x - b.x, a.y - b.y); }
Complex operator*(Complex a, Complex b) { return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }

void fft(int n,Complex * a){
    if(n<=1)
        return ;
    int mid = n >> 1;
    Complex a1[mid], a2[mid];
    for (int i = 0; i < mid;i++){
        a1[i] = a[i << 1];
        a2[i] = a[(i << 1) + 1];
    }
    fft(mid, a1);
    fft(mid, a2);
    Complex w1(cos(PI/mid),sin(PI/mid)),w(1,0),x;
    for (int i = 0; i < mid;i++){
        x = w * a2[i];
        a[i] = a1[i] + x;
        a[i + mid] = a1[i] - x;
        w =w* w1;
    }
}
void ifft(int n,Complex * a){
    if(n<=1)
        return ;
    int mid = n >> 1;
    Complex a1[mid], a2[mid];
    for (int i = 0; i < mid;i++){
        a1[i] = a[i << 1];
        a2[i] = a[(i<<1) + 1];
    }
    ifft(mid, a1);
    ifft(mid, a2);
    Complex w1(cos(PI/mid),-sin(PI/mid)),w(1,0),x;
    for (int i = 0; i < mid;i++){
        x = w * a2[i];
        a[i] = a1[i] + x;
        a[i + mid] = a1[i] -x;
        w =w*w1;
    }
    
}
int main(){
    int n, m;
    cin >> n >> m;
    for (int i = 0; i <=n;i++){
        cin >> a[i].x;
    }
    for (int i = 0; i <=m;i++){
        cin >> b[i].x;
    }
    int k = 1;
    while(k<=n+m)
        k <<= 1;
    fft(k, a);
    fft(k, b);
    for (int i = 0; i <= k; i++){
        a[i] = a[i] * b[i];
    }
    ifft(k, a);
    for (int i = 0; i <=(n+m);i++){
        cout << int(a[i].x / k + 0.5) <<" ";
    }
    //system("pause");
    return 0;
}

给出我的题解，其实就是根据前面给的伪代码写就行。

大整数乘法
如何利用FFT来计算大整数的乘法呢？其实，完全可以把一个多项式看成一个长整数，当 $x = 10$ 带入多项式，多项式的系数就是长整数的数位，那么就完全可以利用多项式的乘法，利用FFT来计算求解。这里贴出一张图可以很快理解其过程，图片内容引用自—大大
请添加图片描述

这里给出一个Leetcode例题：传送门，利用FFT来求大整数的乘法。下面是该题的题解。

class Solution {
    struct Complex{
	    double x, y;
	    Complex(){x=y=0;}
	    Complex(double x,double y){
	        this->x = x;
	        this->y = y;
	    }
	}a[1000], b[1000];
public:
    const double PI = acos(-1);
friend Complex operator+(Complex a, Complex b) { return Complex(a.x + b.x, a.y + b.y); }
friend Complex operator-(Complex a, Complex b) { return Complex(a.x - b.x, a.y - b.y); }
friend Complex operator*(Complex a, Complex b) { return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
    void fft(int n,Complex * a){
    if(n<=1)
        return ;
    int mid = n >> 1;
    Complex a1[mid], a2[mid];
    for (int i = 0; i < mid;i++){
        a1[i] = a[i << 1];
        a2[i] = a[(i << 1) + 1];
    }
    fft(mid, a1);
    fft(mid, a2);
    Complex w1(cos(PI/mid),sin(PI/mid)),w(1,0),x;
    for (int i = 0; i < mid;i++){
        x = w * a2[i];
        a[i] = a1[i] + x;
        a[i + mid] = a1[i] - x;
        w =w* w1;
    }
}
    void ifft(int n,Complex * a){
    if(n<=1)
        return ;
    int mid = n >> 1;
    Complex a1[mid], a2[mid];
    for (int i = 0; i < mid;i++){
        a1[i] = a[i << 1];
        a2[i] = a[(i<<1) + 1];
    }
    ifft(mid, a1);
    ifft(mid, a2);
    Complex w1(cos(PI/mid),-sin(PI/mid)),w(1,0),x;
    for (int i = 0; i < mid;i++){
        x = w * a2[i];
        a[i] = a1[i] + x;
        a[i + mid] = a1[i] -x;
        w =w*w1;
    } 
    }
    string multiply(string num1, string num2) {
        if(num1=="0"||num2=="0") return "0";
        int n =num1.size();int m = num2.size();
        for(int i=0;i<n;i++){
            a[i].x=num1[i]-'0';
        }
        for(int i=0;i<m;i++){
            b[i].x=num2[i]-'0';
        }
        int k=1;
        while(k<=(n+m)) k<<=1;
        fft(k,a);
        fft(k,b);
        for(int i=0;i<=k;i++){
            a[i]=a[i]*b[i];
        }
        ifft(k,a);
        //处理进位
        string ans="";
        int flag=0;
        vector<int>as;
        n--,m--;
        for(int i=(n+m);i>=0;i--){
        as.push_back(int(a[i].x/k+0.5));
        }
        for(auto x :as){
            ans+=to_string((x+flag)%10);
            flag=(x+flag)/10;
        }
        if(flag){
            ans+=to_string(flag);
        }
        //去掉前导0
        /*for(int i=0;i<ans.size();i++){
            if(ans[i]=='0'){
                ans=ans.substr(0,i)+ans.substr(i+1);
                i--;
            }
            else break;
        }*/
        reverse(ans.begin(),ans.end());
        
        return ans;
    }
};

其实直接使用上面求多项式的FFT模板就可以通过，不过在写这题的时候还是遇到了不少坑，要注意细节！

其实归于本质，FFT的加速本质是分治算法的利用。可以参考卜东波老师的讲解，很容易懂，让人醍醐灌顶，录播课在b站上搜得到。

其实，在之前利用DFT求多项式特定点值的时候，本质是利用傅里叶变换转换到频域，然后利用频域的性质，可以直接相乘，而不需要像时域下进行卷积。傅里叶变换能够转换到频域的本质是：多项式系数向量与周期向量做内积。这里的周期向量就是上文中提到的特定点。简单的理解就是，点乘的过程中能够过滤掉其他频率，而求出来的 $y (x)$ 就能够表现出频域的变化。这是利用了 $c o s 和 s i n$ 相乘做积分，周期不同等于0的性质。

说点题外话，卜算子的课确实讲的很好，有一些我之前不理解的点就突然醍醐灌顶了,要是我本科就能听到他的算法课就好了。同时学习的时候回忆起了大一ACM校赛的最后一题，就是利用FFT计算大整数的乘法，同时想起了那年ACM省赛的一道FFT模板题，想起了被通信原理和数字信号处理虐的死去活来的日子，真是感慨万千。

gokingd

关注

6
点赞
踩
11

收藏

觉得还不错? 一键收藏
4
评论
基于FFT的大整数乘法

多项式求值对于多项式f(x)=a0+a1∗x1+a2∗x2+...+an−1∗xn−1f(x)=a_{0}+a_{1}*x^{1}+a_{2}*x^{2}+...+a_{n-1}*x^{n-1}f(x)=a0+a1∗x1+a2∗x2+...+an−1∗xn−1通常我们要花费O(n2)O(n^2)O(n2)的时间来求其nnn个点的处的取值。1965年Cooley and Tukey提出的FFT能够在O(nlogn)O(nlogn)O(nlogn)的时间求得nnn个点处的值。其本质是利用复数的周期性和
复制链接

扫一扫

专栏目录