好久没写算法了，这才是本呀，还得拾起来.
快速傅立叶变换实现两个多项式相乘，求乘积的系数
例如，求(n^2 + 2*n + 1)(2*n^2 + 1)，最高次幂为2
输入： 2 2（两个多项式的最高次幂）
1 2 1（第一个多项式各项系数）
2 0 1（第二个多项式各项系数）
输出：2 4 3 2 1
即：2*n^4 + 4*n^3 + 3*n^2 + 2*n + 1
1:多项式表达法
次数表达法：
次数界n（最高次数为n-1）的多项式: 的系数表达为一个由系数组成的向量a=（a0,a1,a2,….an-1）
点值表达法:
把多项式理解成一个函数，然后用函数上的点表示这个函数
（两点确定一条直线，三点确定一个二次函数…n+1个点确定一个n次函数，以此类推）
次数界为n的多项式的点值表达法就是一个n个点对组成的集合
2:单位复数根
如果一个数的n次方能够变回1，那么我们叫他n次复数根，记作w_n^k
n次单位复数根刚好有n个, 他们是e^(2kπ/n i), k=0,1,2…n−1, 关于复数指数的定义如下：
e^ui=cos〖(u)+sin(u)〗 i。
他们均匀的分布在了这个圆上，离原点的距离都是1
下图以四次单位复数根为例 
关于复数根的几个定理和引理：
消去引理： 对任何整数 n≥0,k≥0,d≥0n≥0,k≥0,d≥0有ωdkdn=ωknωdndk=ωnk
证明： ωdkdn=(e2iπdn)dk=(e2iπn)k=ωknωdndk=(e2iπdn)dk=(e2iπn)k=ωnk
一个推论： 对任意偶数 n > 0 有 ωn2n=ω2=−1ωnn2=ω2=−1
证明：设n = 2*k那么 ωn2n=ωk2k=ω2=eπ=cos(π)+sin(π)i=−1ωnn2=ω2kk=ω2=eπ=cos(π)+sin(π)i=−1
折半引理：如果n > 0是偶数, 那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合
证明：实际上这个引理就是证明了(ωk+n2n)2=ω2k+nn=ω2knωnn=(ωkn)2(ωnk+n2)2=ωn2k+n=ωn2kωnn=(ωnk)2
折半引理对于采用分治对多项式系数向点值表达的转换有很大作用, 保证了递归的子问题是原问题规模的一半
求和引理：对任意整数n≥1n≥1和不能被n整除的非负整数k, 有
Σj=0n−1(ωkn)j=0Σj=0n−1(ωnk)j=0
这个问题通过等比数列求和公式就可以得到：Σn−1j=0(ωkn)j=(ωkn)k−1ωkn−1=1k−1ωkn−1=0
DFT
在DFT变换中, 希望计算多项式A(x)在复数根ω0n,ω1n,...,ωn−1nωn0,ωn1,...,ωnn−1处的值, 也就是求
yk=A(ωkn)=Σj=0n−1ajωkjnyk=A(ωnk)=Σj=0n−1ajωnkj
称向量y=(y0,y1,...,yn−1)y=(y0,y1,...,yn−1)是系数向量a=(a0,a1,...,an−1)a=(a0,a1,...,an−1)的离散傅里叶变换, 记为y=DFTn(a)y=DFTn(a)
FFT
直接计算DFT的复杂度是O(n2)O(n2), 而利用复数根的特殊性质的话, 可以在O(nlogn)O(nlogn)的时间内完成, 这个方法就是FFT方法, 在FFT方法中采用分治策略来进行操作, 主要利用了消去引理之后的那个推论
在FFT的策略中, 多项式的次数是2的整数次幂, 不足的话再前面补0, 每一步将当前的多项式A(x), 次数是2的倍数, 分成两个部分：
A[0](x)=a0+a2x+a4x2+...+an−2xn2−1A[0](x)=a0+a2x+a4x2+...+an−2xn2−1
A[1](x)=a1+a3x1+a5x2+...+an−1xn2−1A[1](x)=a1+a3x1+a5x2+...+an−1xn2−1
于是就有了
A(x)=A[0](x2)+xA[1](x2)A(x)=A[0](x2)+xA[1](x2)
那么我们如果能求出次数界是n2n2的多项式A[0](x)A[0](x)和A[1](x)A[1](x)在n个n次单位复数根的平方处的取值就可以了, 即在
(ω0n)2,(ω1n)2,(ω2n)2,...,(ωn−1n)2(ωn0)2,(ωn1)2,(ωn2)2,...,(ωnn−1)2
处的值, 那么根据折半引理, 这n个数其实只有n2n2个不同的值, 也就是说, 对于每次分出的两个次数界n2n2的多项式, 只需要求出其n2n2个不同的值即可, 那么问题就递归到了原来规模的一半, 也就是说如果知道了两个子问题的结果, 当前问题可以在两个子问题次数之和的复杂度内解决, 那么这样递归问题的复杂度将会是O(nlogn)
#include <iostream>
using namespace std;
#include <complex>//complex是复数的意思，C++中已经提供了实现复数的类complex，而complex<double>则是声明一些实部和虚部都为double类型的复数变量。
#include <cmath>
double PI = acos(-1);
complex<double> a[400010], b[400010], c[400010];
void fft(complex<double> *a, int n, int op)complex是复数的意思，C++中已经提供了实现复数的类complex，而complex<double>则是声明一些实部和虚部都为double类型的复数变量。
#include <cmath>
double PI = acos(-1);
complex<double> a[400010], b[400010], c[400010];
void fft(complex<double> *a, int n, int op){
if (n == 1) return;
complex<double> w(1, 0), wn(cos(2*PI/n), sin(2*PI*op/n)), a1[n>>1], a2[n>>1];
for (int i = 0; i < (n>>1); i++)
a1[i] = a[i<<1], a2[i] = a[(i<<1)+1];
fft(a1, n>>1, op), fft(a2, n>>1, op);
for (int i = 0; i < (n>>1); i++, w*=wn)
a[i] = a1[i] + w*a2[i], a[i+(n>>1)] = a1[i] - w*a2[i];
}
int main()
{
int n, m;
scanf("%d%d", &n, &m);//两个式子中最高项的次数
for (int i = 0; i <= n; i++) scanf("%lf", &a[i]);//依次输入每一项的系数
for (int i = 0; i <= m; i++) scanf("%lf", &b[i]);
m += n, n = 1;
while (n <= m) n <<= 1;
fft(a, n, 1), fft(b, n, 1);
for (int i = 0; i < n; i++) c[i] = a[i] * b[i];
fft (c, n, -1);
for (int i = 0; i <= m; i++) printf("%.0lf ", floor(c[i].real()/n));
return 0;
}
<img data-cke-saved-src="https://img-blog.csdnimg.cn/2022010613412770032.png" src="https://img-blog.csdnimg.cn/2022010613412770032.png" alt="" />
FFT递归版：
#include <iostream>
using namespace std;
#include <complex>
#include <cmath>
double PI = acos(-1);
complex<double> a[400010], b[400010], c[400010];
void fft(complex<double> *a, int n, int op)
{
if (n == 1) return;
complex<double> w(1, 0), wn(cos(2*PI*op/n), sin(2*PI*op/n)), a1[n>>1], a2[n>>1];
for (int i = 0; i < (n>>1); i++)
a1[i] = a[i<<1], a2[i] = a[(i<<1)+1];
fft(a1, n>>1, op), fft(a2, n>>1, op);
for (int i = 0; i < (n>>1); i++, w*=wn)
a[i] = a1[i] + w*a2[i], a[i+(n>>1)] = a1[i] - w*a2[i];
}
int main()
{
int n, m;
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; i++) scanf("%lf", &a[i]);
for (int i = 0; i <= m; i++) scanf("%lf", &b[i]);
m += n, n = 1;
while (n <= m) n <<= 1;
fft(a, n, 1), fft(b, n, 1);
for (int i = 0; i < n; i++) c[i] = a[i] * b[i];
fft (c, n, -1);
for (int i = 0; i <= m; i++) printf("%d ", int(c[i].real()/n+0.5));
return 0;
}<img data-cke-saved-src="https://img-blog.csdnimg.cn/2022010613412852465.png" src="https://img-blog.csdnimg.cn/2022010613412852465.png" alt="" />
FFT迭代版：
#include <iostream>
using namespace std;
#include <cmath>
double PI = acos(-1);
const int MAXN = 4*1e5+10;
int r[MAXN];
struct Complex{
double r, i;
Complex(){}
Complex(double _r, double _i){ r = _r, i = _i; }
Complex operator +(const Complex &y) { return Complex(r+y.r, i+y.i); }
Complex operator -(const Complex &y) { return Complex(r-y.r, i-y.i); }
Complex operator *(const Complex &y) { return Complex(r*y.r - i*y.i, r*y.i+i*y.r); }
Complex operator *=(const Complex &y) { double t = r; r = r*y.r - i*y.i, i = t*y.i + i*y.r; }
}a[MAXN], b[MAXN];
void fft(Complex *a, int len, int op)
{
Complex tt;
for (int i = 0; i < len; i++) if (i < r[i])
tt = a[i], a[i] = a[r[i]], a[r[i]] = tt;
for (int i = 1; i < len; i <<= 1)
{
Complex wn(cos(PI/i), sin(PI*op/i));
for (int k=0, t=(i<<1); k < len; k += t)
{
Complex w(1, 0);
for (int j = 0; j < i; j++, w *= wn)
{
Complex t = w*a[k+j+i];
Complex u = a[k+j];
a[k+j] = u + t;
a[k+j+i] = u - t;
}
}
}
if (op == -1) for (int i = 0; i < len; i++) a[i].r /= len, a[i].i /= len;
}
int main()
{
int n, m, L, i, x;
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; i++) scanf("%lf", &a[i]);
for (int i = 0; i <= m; i++) scanf("%lf", &b[i]);
m += n;
for (n=1, L=0; n <= m; n <<= 1) L++;
for (i=0, x=L-1; i < n; i++) r[i] = (r[i>>1]>>1) | ((i&1)<<x);
fft(a, n, 1), fft(b, n, 1);
for (int i = 0; i < n; i++) a[i] *= b[i];
fft (a, n, -1);
for (int i = 0; i <= m; i++) printf("%d ", int(a[i].r+0.5));
return 0;
}<img data-cke-saved-src="https://img-blog.csdnimg.cn/2022010613412828763.png" src="https://img-blog.csdnimg.cn/2022010613412828763.png" alt="" />
总结：普通的
扫一扫
1835
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
