快速傅里叶变换FFT
大整数乘法计算的发展历程是缓慢的...在朴素的大整数乘法计算和多项式乘法中,两个多项式f(x), g(x)的任意两项之间都要相乘一次,时间复杂度为O(n^2)
。后来出现了用分治法计算大整数乘法,把时间复杂度优化到O(n^{1.585})
,然后经历了相当长的时间后才发明出了更快的计算方法--快速傅里叶变换(Fast Fourier Transform),时间复杂度为O(nlogn)
,为目前计算大数乘法最快的算法。
多项式的表示
一个多项式有两种表示方法,系数表示法和点值表示法。系数表示法就是我们最常用的用来表示一个多项式的方法,即多项式每个项的系数组成的集合。一个最高n次的多项式可以用n+1个系数来表示,也可以用n+1个横坐标不同的点来表示。n+1个不同的点可以唯一表示一个最高n次多项式。
给定两个多项式f(x), g(x),现计算其乘法h(x) = f(x)*g(x)
,若这两个多项式是用系数表示法表示的,那么需要O(n^2)
的时间去计算;若是用点值表示法表示的,假设取到的点的横坐标的都是相同的,即(x1, f(x1)), (x2, f(x2)), (x3, f(x3))....(xn, f(xn))和(x1, g(x1)), (x2, g(x2)), (x3, g(x3))....(xn, g(xn)),那么仅需要O(n)的时间即可计算出h(x)的点值表示法(x1, f(x1)*g(x1)), (x2, f(x2)*g(x2)), (x3, f(x3)*g(x3))....(xn, f(xn)*g(xn))。可以利用这一过程来加速计算多项式乘法。
难点是如何将多项式从系数表示法转化为点值表示法。我们知道计算一个n次多项式的一个f(x)需要O(n)的时间,那转换成点值表示法岂不是也需要O(n^2)
的时间。。快速傅里叶变换做的工作就是快速进行转换这个工作。
这里需要引入复数。
复数
如何快速从系数表示法转化为点值表示法呢?首先需要选择n+1个横坐标,这里选择x^n=1
在复数意义上的n个单位根来表示,利用单位根的一些性质来快速完成转化!
一个复数a+bi
可以像向量一样在二维坐标里表示为坐标(a, b),也可以像极坐标一样用一对(模角,幅角)来表示,模角是复数的长度,幅角是该复数和x正半轴的夹角。两复数相乘时满足:幅角相乘,模角相加。所以x^n=1
的n个根其实就是把一个圆心在原点的单位元n等分,得到的n个点。
设w_{n}^{k}
表示n次复根的第k个单位根。w_{n}^{k}
对应的坐标为(cos(\frac{k}{n}2\pi), sin(\frac{k}{n}2\pi))
,也就是复数cos(\frac{k}{n}2\pi)+i*sin(\frac{k}{n}2\pi)
,根据欧拉公式也可表示成e^{i·\frac{k}{n}2\pi}
。用欧拉公式可以很容易的解释为什么两个复数相乘结果是模角相乘,幅角相加。
三个重要性质:
w_{n}^{n} = w_{n}^{0} = 1
w_{n}^{k} = w_{2n}^{2k}
w_{n}^{k+\frac{n}{2}} = -w_{n}^{k}
这三个性质都很容易推出来。
第三个性质中,因为w_{n}^{\frac{n}{2}}
刚好是-1,所以去掉之后还需要再乘以-1.
蝴蝶变换
如何把复数应用进去应该是FFT最难理解的部分了。
我们现在要做的是:对于一个n-1次多项式f(x),把n个单位复根带入,求得n个复值,如何快速的求出。
我们把f(x)的奇次项的系数和偶此项的系数分别表示成一个多项式,即
A_0(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{\frac{n}{2}-1}
A_1(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{\frac{n}{2}-1}
那么就有f(x) = A_0(x^2) + x * A_1(x^2)
.
我们再将n次单位复根带入,设k<=\frac{n}{2}
,得
f(w_{n}^{k}) = A_0(w_{n}^{2k}) + w_{n}^{k}*A_1(w_{n}^{2k}) = A_0(w_{\frac{n}{2}}^{k}) + w_{n}^{k}*A_1(w_{\frac{n}{2}}^{k})
f(w_{n}^{k+\frac{n}{2}}) = A_0(w_{n}^{2k+n}) + w_{n}^{k+\frac{n}{2}}*A_1(w_{n}^{2k+n}) = A_0(w_{\frac{n}{2}}^{k}) - w_{n}^{k}*A_1(w_{\frac{n}{2}}^{k})
f(x)是n-1次多项式,需要计算n次单位复根带入得到的值;A_0(x), A_1(x)
是\frac{2}{n}-1
次多项式,已经计算出将\frac{2}{n}
次单位复根带入得到的值,这样就可以在A_0(x), A_1(x)
的基础上线性求出f(x)的n个复根!A_0(x), A_1(x)
也可以用类似的求法递归求下去,共需要求logn
层,所以总的时间复杂度是O(nlogn)
.
注意到每层循环中n都必须可以被2整除,所以第一层的n必须是2的幂,若所给的多项式的最高次不是2的幂,那么就补0.
快速傅里叶逆变换IFFT
得到目标多项式的点值表示法之后,我们还要把它转化成系数表示法,这一过程叫快速傅里叶逆变换。
把f(x)转化成点值表示法后,再把这些点值作为另一个多项式的系数,并把单位复根的倒数作为根代入,得到的值就是f(x)的系数,所以这里相当于是再做一次FFT即可。
代码:
const double PI = acos(-1.0);
void change(Complex y[], int len){
for (int i=1, j=len/2, k; i<len-1; i++){
if (i<j) swap(y[i], y[j]);
k = len/2;
while (j>=k){
j -= k;
k /= 2;
}
if (j<k) j += k;
}
}
void FFT(Complex y[], int len, int on){ //快速傅里叶变换,y表示要带入的len个复数,on表示是求正变换还是逆变换。
change(y, len); //直接化成递归的最底一层,即最低一层中各项系数的位置。
for (int h=2; h<=len; h<<=1){ //h表示该层每个函数有h个元素,开始向上迭代
Complex wn(cos(PI*2/h), sin(PI*2*on/h)); //wn为第1个单位根/单位根倒数
for (int j=0; j<len; j+=h){ //每层的第一个元素
Complex w(1, 0); //w表示第几单位根,求完一次后,都变成下一个单位根。
for (int k=j; k<j+h/2; k++){
Complex u = y[k], t = w*y[k+h/2];
y[k] = u+t, y[k+h/2] = u-t;
w = w*wn;
}
}
}
if (on==-1){ //若为逆变换,还需要再除以n才能得到原多项式的系数。
for (int i=0; i<len; i++) y[i].x /= len;
}
}
复数的使用:
struct Complex{
double x, y; //实部和虚部
Complex(double a = 0.0, double b = 0.0){
x = a, y = b;
}
Complex operator + (const Complex &a) const {
return Complex(x+a.x, y+a.y);
}
Complex operator - (const Complex &a) const {
return Complex(x-a.x, y-a.y);
}
Complex operator * (const Complex &a) const { //简单推出
return Complex(x*a.x-y*a.y, x*a.y+y*a.x);
}
};
题目
P1919 【模板】A*B Problem升级版
大数相乘不同于多项式乘法的地方是在计算完成后需要进位。
#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e6+500000;
const double PI = acos(-1.0);
struct Complex{
double x, y;
Complex(double a = 0.0, double b = 0.0){
x = a, y = b;
}
Complex operator + (const Complex &a) const {
return Complex(x+a.x, y+a.y);
}
Complex operator - (const Complex &a) const {
return Complex(x-a.x, y-a.y);
}
Complex operator * (const Complex &a) const {
return Complex(x*a.x-y*a.y, x*a.y+y*a.x);
}
};
char s1[maxn], s2[maxn];
struct Complex x1[maxn*2], x2[maxn*2];
int len;
int ans[maxn*2];
void change(Complex y[], int len){
for (int i=1, j=len/2, k; i<len-1; i++){
if (i<j) swap(y[i], y[j]);
k = len/2;
while (j>=k){
j -= k;
k /= 2;
}
if (j<k) j += k;
}
}
void FFT(Complex y[], int len, int on){
change(y, len); //直接化成递归的最底一层
for (int h=2; h<=len; h<<=1){ //h表示该层每个函数有h个元素
Complex wn(cos(PI*2/h), sin(PI*2*on/h)); //
for (int j=0; j<len; j+=h){ //每层的第一个元素
Complex w(1, 0); //w表示第几单位根
for (int k=j; k<j+h/2; k++){
Complex u = y[k], t = w*y[k+h/2];
y[k] = u+t, y[k+h/2] = u-t;
w = w*wn;
}
}
}
if (on==-1){
for (int i=0; i<len; i++) y[i].x /= len;
}
}
int main(){
while (~scanf ("%s %s", s1, s2)){
int len1 = strlen(s1), len2 = strlen(s2);
len = 1;
while (len < len1+len2) len <<= 1;
for (int i=0; i<len1; i++) x1[i] = Complex(s1[len1-i-1]-'0', 0.0);
for (int i=0; i<len2; i++) x2[i] = Complex(s2[len2-i-1]-'0', 0.0);
for (int i=len1; i<len; i++) x1[i] = Complex(0, 0);
for (int i=len2; i<len; i++) x2[i] = Complex(0, 0);
FFT(x1, len, 1);
FFT(x2, len, 1);
for (int i=0; i<len; i++) x1[i] = x1[i] * x2[i];
FFT(x1, len, -1);
for (int i=0; i<len; i++) ans[i] = int(x1[i].x + 0.5);
for (int i=0; i<len; i++){
ans[i+1] += ans[i] / 10;
ans[i] %= 10;
}
while (ans[len-1]==0 && len>1) len--;
for (int i=len-1; i>=0; i--) putchar(ans[i]+'0');
putchar('\n');
}
return 0;
}
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 1e6+500000;
const double PI = acos(-1.0);
struct Complex{
double x, y;
Complex(double a = 0.0, double b = 0.0){
x = a, y = b;
}
Complex operator + (const Complex &a) const {
return Complex(x+a.x, y+a.y);
}
Complex operator - (const Complex &a) const {
return Complex(x-a.x, y-a.y);
}
Complex operator * (const Complex &a) const {
return Complex(x*a.x-y*a.y, x*a.y+y*a.x);
}
};
struct Complex x1[maxn*2], x2[maxn*2];
int len;
ll ans[maxn*2];
void change(Complex y[], int len);
void FFT(Complex y[], int len, int on);
int main(){
int len1, len2, tmp;
while (~scanf ("%d %d", &len1, &len2)){
len = 1, len1++, len2++;
while (len < len1+len2) len <<= 1;
for (int i=0; i<len1; i++){
scanf ("%d", &tmp);
x1[i] = Complex(tmp, 0.0);
}
for (int i=0; i<len2; i++){
scanf ("%d", &tmp);
x2[i] = Complex(tmp, 0.0);
}
for (int i=len1; i<len; i++) x1[i] = Complex(0, 0);
for (int i=len2; i<len; i++) x2[i] = Complex(0, 0);
FFT(x1, len, 1);
FFT(x2, len, 1);
for (int i=0; i<len; i++) x1[i] = x1[i] * x2[i];
FFT(x1, len, -1);
for (int i=0; i<len; i++) ans[i] = ll(x1[i].x + 0.5);
while (ans[len-1]==0 && len>1) len--;
printf ("%lld", ans[0]);
len1--, len2--;
for (int i=1; i<=len1+len2; i++) printf (" %lld", ans[i]);
}
return 0;
}
void change(Complex y[], int len){
for (int i=1, j=len/2, k; i<len-1; i++){
if (i<j) swap(y[i], y[j]);
k = len/2;
while (j>=k){
j -= k;
k /= 2;
}
if (j<k) j += k;
}
}
void FFT(Complex y[], int len, int on){
change(y, len); //直接化成递归的最底一层
for (int h=2; h<=len; h<<=1){ //h表示该层每个函数有h个元素
Complex wn(cos(PI*2/h), sin(PI*2*on/h)); //
for (int j=0; j<len; j+=h){ //每层的第一个元素
Complex w(1, 0); //w表示第几单位根
for (int k=j; k<j+h/2; k++){
Complex u = y[k], t = w*y[k+h/2];
y[k] = u+t, y[k+h/2] = u-t;
w = w*wn;
}
}
}
if (on==-1){
for (int i=0; i<len; i++) y[i].x /= len;
}
}
0 条评论