多项式乘法 快速傅里叶变换

好久没写算法了，这才是本呀，还得拾起来.

快速傅立叶变换实现两个多项式相乘，求乘积的系数

例如，求(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 有 

证明：设n = 2*k那么 

折半引理：如果n > 0是偶数, 那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合

证明：实际上这个引理就是证明了

折半引理对于采用分治对多项式系数向点值表达的转换有很大作用, 保证了递归的子问题是原问题规模的一半

求和引理：对任意整数和不能被n整除的非负整数k, 有

这个问题通过等比数列求和公式就可以得到：

DFT

在DFT变换中, 希望计算多项式A(x)在复数根处的值, 也就是求

称向量是系数向量的离散傅里叶变换, 记为

FFT

直接计算DFT的复杂度是, 而利用复数根的特殊性质的话, 可以在的时间内完成, 这个方法就是FFT方法, 在FFT方法中采用分治策略来进行操作, 主要利用了消去引理之后的那个推论

在FFT的策略中, 多项式的次数是2的整数次幂, 不足的话再前面补0, 每一步将当前的多项式A(x), 次数是2的倍数, 分成两个部分：

于是就有了

那么我们如果能求出次数界是的多项式和在n个n次单位复数根的平方处的取值就可以了, 即在

处的值, 那么根据折半引理, 这n个数其实只有个不同的值, 也就是说, 对于每次分出的两个次数界的多项式, 只需要求出其个不同的值即可, 那么问题就递归到了原来规模的一半, 也就是说如果知道了两个子问题的结果, 当前问题可以在两个子问题次数之和的复杂度内解决, 那么这样递归问题的复杂度将会是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)

{
	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;
}

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;
}

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;
}

总结：普通的傅里叶变换用时19.15秒，FFT递归版用时19.77秒，FFT迭代版用时7.005秒，可见不同方法的处理对程序运行时间还是 有差别的。