基于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个点处的值。其本质是利用复数的周期性和分治算法消除计算的冗余。故FFT加速求值时，只能计算特殊点处的取值，如 1,w,w2,...,wn−11,w,w^2,...,w^{n-1}1,w,w2,...,wn−1,其中w=e2π∗i/nw=e^{2\pi*i/n}w=e2π∗i/n
具体过程这里不做推导。直接给出FFT计算多项式值的伪代码。

–引用卜算的PPT

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

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

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

多项式乘法
对于多项式A(x)A(x)A(x)和B(x)B(x)B(x)求C(x)C(x)C(x)
A(x)=a0+a1∗x1+a2∗x2+...+an−1∗xn−1B(x)=b0+b1∗x1+b2∗x2+...+bn−1∗xn−1C(x)=A(x)∗B(x)=c0+c1∗c1+c2∗x2+...+c2n−2∗x2n−2A(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} A(x)=a0+a1∗x1+a2∗x2+...+an−1∗xn−1B(x)=b0+b1∗x1+b2∗x2+...+bn−1∗xn−1C(x)=A(x)∗B(x)=c0+c1∗c1+c2∗x2+...+c2n−2∗x2n−2
ck=Σi=0kaibk−ic_{k}=\Sigma_{i=0}^{k}a_{i}b_{k-i}ck=Σi=0kaibk−i的每一项其实是计算卷积，求出C需要的时间为O(n2)O(n^2)O(n2)。这时候就要放上一张经典图了。

从我们上面的逻辑来看，就是先通过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=10x=10x=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)y(x)y(x)就能够表现出频域的变化。这是利用了cos和sincos和sincos和sin相乘做积分，周期不同等于0的性质。

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

基于FFT的大整数乘法相关推荐

使用快速傅里叶变换计算大整数乘法-代码
在上一篇随笔"使用快速傅里叶变换计算大整数乘法"中,已经讲述了使用快速傅里叶变换计算大整数乘法的原理.在这一篇随笔中,我们就使用快速傅里叶变换来实现一个提供任意精度的算术运算的静态 ...
大整数乘法（Karatsuba算法的字符串形式的C++实现）
#include <iostream> #include <sstream> #include <cstring> using namespace std;/函数声 ...
【FFTNTT入门】大整数乘法
问题:给定两个大整数 A A A 和 B B B, A A A 和 B B B 的长度为 n n n 和 m m m,求 A A A 和 B B B 的乘积 1. 朴素做法思考小学数学中两个数的乘法 ...
大整数乘法c语言代码_大整数乘法
大整数乘法和我们小学学过的乘法公式一样(如下图),就是按位相乘,两个数中的每一位彼此相乘,然后将相同列的结果加起来,最后统一处理进位即可. #include <iostream> #inc ...
信息学奥赛一本通 1307：【例1.3】高精度乘法 | 1174：大整数乘法 | OpenJudge NOI 1.13 09:大整数乘法
[题目链接] ybt 1307:[例1.3]高精度乘法 ybt 1174:大整数乘法 OpenJudge NOI 1.13 09:大整数乘法 [题目考点] 1. 高精度考察:高精乘高精高精度计算讲 ...
信息学奥赛一本通（1174：大整数乘法）
1174:大整数乘法时间限制: 1000 ms 内存限制: 65536 KB 提交数: 12480 通过数: 7002 [题目描述] 求两个不超过200位的非负整数的积. [ ...
大整数乘法--leetcode Multiply Strings
大整数乘法本文转载自http://www.cnblogs.com/TenosDoIt/p/3735309.html 我们在日常的大整数计算中,通常是把它转化为字符型计算.这道题的思路就和我们小学计算 ...
python两数相乘代码_Python 实现大整数乘法算法的示例代码
我们平时接触的长乘法,按位相乘,是一种时间复杂度为 O(n ^ 2) 的算法.今天,我们来介绍一种时间复杂度为 O (n ^ log 3) 的大整数乘法(log 表示以 2 为底的对数). 介绍原理 ...
计算机算法设计与分析大整数乘法
大整数乘法问题描述求两个不超过200位的非负整数的积. 输入形式有两行,每行是一个不超过200位的非负整数,没有多余的前导0. 输出形式一行,即相乘后的结果.结果里不能有多余的前导0,即如果结 ...

基于FFT的大整数乘法

基于FFT的大整数乘法相关推荐

最新文章

热门文章