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

2019年11月7日 96次阅读

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

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

例如，求(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≥0 n≥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≥1 n≥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ωkjn yk=A(ωnk)=Σj=0n−1⁡ajω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−1 A[0](x)=a0+a2x+a4x2+…+an−2xn2−1

A[1](x)=a1+a3x1+a5x2+...+an−1xn2−1 A[1](x)=a1+a3x1+a5x2+…+an−1xn2−1

于是就有了

A(x)=A[0](x2)+xA[1](x2) A(x)=A[0](x2)+xA[1](x2)

那么我们如果能求出次数界是
n2 n2
的多项式
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个数其实只有
n2 n2
个不同的值, 也就是说, 对于每次分出的两个次数界
n2 n2
的多项式, 只需要求出其
n2 n2
个不同的值即可, 那么问题就递归到了原来规模的一半, 也就是说如果知道了两个子问题的结果, 当前问题可以在两个子问题次数之和的复杂度内解决, 那么这样递归问题的复杂度将会是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秒，可见不同方法的处理对程序运行时间还是 有差别的。

点赞