算法导论第三十(30)章多项式与快速傅里叶变换

最新推荐文章于 2024-09-29 15:06:28 发布

Flammable_ice

最新推荐文章于 2024-09-29 15:06:28 发布

阅读量9.8k

点赞数 6

分类专栏：《算法导论》

本文链接：https://blog.csdn.net/z84616995z/article/details/24308087

版权

本文介绍了如何使用快速傅里叶变换（FFT）进行多项式乘法，包括点值乘法、系数转换和逆变换的过程。通过实例详细解析了30.1-1至30.1-5的练习，探讨了多项式的确定性、计算复杂度以及笛卡尔和的计算方法。同时，讲解了DFT、FFT的基本概念、引理和应用，包括DFT逆的计算。

摘要由CSDN通过智能技术生成

由于在第三十一章数论算法中遇到几个关于超大数乘法的问题促使我需要学这章，具体请看第三十一章练习31.1-8，31.1-12与31.1-13.

基本概念：

大整数的加减乘除

系数形式表示的多项式的快速乘法

1.两个n次多项式a与b分别计算在2n次单位复数根下的对应多项式值。计算过程使用的是FFT。所以时间为O(nlgn)

2.随后进行点值乘法，根据c(w0)=a(w0)*b(w0)在n次单位复数根下的对应的多项式值c(w0),那么记录点c(w0,c(w0)) 以此类推，记录一系列C=(c0,c1....cn).，其时间为O(n)

3.把一系列C在单位复数根下的点值表达通过DFT逆计算方式转换成多项式C的系数表达，由此求得结果。其时间为O(nlgn),所以总的运行时间是O(nlgn)+O(nlgn)+O(n)=O(nlgn)

30.1-1 运用等式(30,1)和(30.2)把下列两个多项式相乘：A(x)=7x^3-x^2+x-10 B(x)=8x^3-6x+3

由于过于简单,直接上结果了。C(x)=56x^6-8x^5-34x^4-53x^3-9x^2+63x-30.

30.1-2 求一个次数界为n的多项式A(x)在某给定点x0的值存在另外一种方法；把多项式A(x)除以多项式(x-x0),得到一个次数界为n-1的商多项式q(x)和余项r,满足A(x)=q(x)(x-x0)+r.很明显，A(x0)=r.请说明如何根据x0和A的系数，在O(n)的时间复杂度内计算出余项r以及q(x)的系数。

//30.1-2另类的计算多项式在某点值方法
#include <iostream>
#include <vector>
using namespace std;
struct Polynomial
{
	int coef;//系数
	int Index;//指数
};
int Division_with_remainder(vector<Polynomial>a,vector<Polynomial>&qx,int x0,const int n)
{
//	vector<Polynomial>temp;
	int temp=0,i=0,s=1,k=a[n-1].Index;
	Polynomial q;
    if (n==1)//如果就一项，那么直接求值即可
    {
		i=a[0].Index;
		q.coef=-a[0].coef;
		q.Index=a[0].Index;
		qx.push_back(q);
		while (i)
		{
			s*=x0;
			i--;
		}
		return s*a[0].coef;
    }
	else//如果多余2项，那么用类似小学的除法式子去除
	{
		while(k>=1)
		{
			q.coef=-a[k].coef*x0;
			cout<<q.coef<<endl;
            a[k-1].coef=a[k-1].coef-a[k].coef*x0;
			q.Index=a[k-1].Index;
			qx.push_back(q);
			k--;
		}
		return a[k].coef;
	}
}
void main()
{
	cout<<"请输入需要在哪个点求值："<<endl;
	int x0,i=0;
	cin>>x0;//由于程序并不非常强壮，所以请不要输入字符或者直接输入ctrl+z,只能输入一些整数。输入系数时也是一样的。
	cout<<"请按由低次到高次依次输入每项系数，且每项次数均不同(注意系数为0的项不能省略输入，直接补0)："<<endl;
	vector<Polynomial>a,qx;//a表示多项式系数向量。
	struct Polynomial x;
	while (cin>>x.coef)//输入多项式系数，从低次项到高次项依次输入，输入ctrl+z结束输入。
	{
		x.Index=i;
		a.push_back(x);
		i++;
	}
	const int n=i;
	cout<<"您输入的多项式为："<<endl;
	while (i)
	{
		if(i==1)cout<<a[i-1].coef<<"x^"<<a[i-1].Index;
		else cout<<a[i-1].coef<<"x^"<<a[i-1].Index<<"+";
		i--;
	}
	cout<<endl;
	cout<<"多项式在点"<<x0<<"的值为："<<Division_with_remainder(a,qx,-x0,n)<<endl;
	int k=0;
	cout<<"多项式的商q(x)=";
	while (k!=n-1)
	{
		if(k==n-2)cout<<qx[k].coef<<"x^"<<qx[k].Index<<endl;
		else cout<<qx[k].coef<<"x^"<<qx[k].Index<<"+";
		k++;
	}
}

30.1-3从A(x)=∑ajx^j的点值表达推导出A^res(x)=∑a(n-1-j)x^j(res是将多项式系数逆转)的点值表达，假设没有一个点是0。

这里简要的说明下思路：先用DFT逆求出A(x)在单位根下的系数表达，然后逆转系数(倒序)再次用FFT求出逆转多项式在单位根下的点值表达。运行时间为O(nlgn).

30.1-4证明：为了唯一确定一个次数界为n的多项式，n个相互不同的点值对是必须的。也就是说，如果给定少于n个不同点值对，则他们不能唯一确定一个次数界为n的多项式。(提示：利用定理30.1，加入一个任意选择的点值对到一个已有n-1个点值对的集合，看看会发生什么？)

如果少于n个不同点值对，那么这n个点里必有相同的点使：xi=xj,那么范德蒙德矩阵的行列式=0=>矩阵是不可逆的=>

所以不能通过矩阵方程确定唯一的系数向量a，得出结论。

30.1-5 说明如何利用等式（30.5）在O(n^2)时间复杂度内进行插值运算。(提示：首先计算多项式π(x-xj)的系数表达，然后把每个项除以(x-xk);参见练习30.1.2.你可以在O(n)时间内计算n个分母的每一个。)

这里就说下思路：假设已知n个点(x0,y0),(x1,y1)...(xn,yn).

1.将n个单位复数根依次带入π(x-xj)得出n个多项式值yi,这个过程为O(n^2).

2.然后利用DFT逆求出它的系数表达，这个过程需要O(nlgn).

3.然后再写一个求yk/π(xj-xk)的循环，这个循环大约O(n^2),把每个求得的值存放在数组a中.

4.于是利用练习30.1-2，对这个求出系数的多项式依次除以(x-xj)并且乘以a[i]，每次做多项式除法需要O(n)时间，然后同时循环n次求累加和，最终这个双重循环需要O(n^2).

最后计算总的运行时间为T(n)=O(n^2)+O(nlgn)+O(n^2)+O(n^2)=O(n^2).

30.1-6 此题没有思路。网上也没找到答案。这道题主要讨论的是多项式的(非)确定性。

30.1-7考虑两个集合A和B，每个集合包含取值范围在0~10n之间的n个整数。我们希望计算出A与B的笛卡尔和，定义如下：C={x+y:x∈A，y∈B}注意到，C中整数值的范围在0~20n之间。我们希望找到C中的元素，并且求出C中的每个元素可表示为A中元素与B中元素和的次数。请在O(nlgn)时间内解决问题.(提示：请用次数至多是10n的多项式来表示A和B)。

多项式A与B是系数均为1的，指数为0~10n之间的n个整数.计算C(x)=A(x)*B(x)可以用DFT在O(nlgn)时间内计算出。然后再遍历一下C，C的系数就是A与B的笛卡尔和项的出现次数

30.2 DFT与FFT

单位复数根：n次单位复数根是满足w^n=1的复数w. 复数的指数形式定义:e^iu=cos(u)+isin(u).

n次单位复数根：wn=e^(2πi/n)。

30.3消去引理：对任何整数n≥0，k≥0，以及d>0, w(dn,dk)=w(n,k)//w的第一个参数是下标，第二个参数是次数。

推论30.4 对任意偶数n>0,有w(n,n/2)=w(2,1)=-1.

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

引理30.6求和引理：对任意整数n≥1和不能被n整除的非负整数k,有∑(w(n,k)^j)=0

DFT:就是求在单位复数根处的多项式的值。

FFT：就是在O(nlgn)时间内计算出DFT(a).

DFT逆：就是把书上伪代码a与y互换，wn^-1替换wn，并将结果÷n即可求得答案。

书上伪代码翻译成C++程序如下：

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <vector>
using namespace std;
#define pi 3.1415926535 
int m=0;
struct  Complex
{
	double real;
	double imag;
	int num;
	Complex(int i=0,int j=0){real=i,imag=j;}
};
Complex q(1,0);
Complex operator+ (Complex const & lhs, Complex const & rhs)
{
	Complex temp;
	temp.real =lhs.real+rhs.real;
	temp.imag = lhs.imag+ rhs.imag;
	return temp;
}
Complex operator- (Complex const & lhs, Complex const & rhs)
{
	Complex temp;
	temp.real = lhs.real - rhs.real;
	temp.imag = lhs.imag - rhs.imag;
	return temp;
}
Complex operator* (Complex const & lhs, Complex const & rhs)
{
	Complex temp;
	temp.real = lhs.real * rhs.real - lhs.imag * rhs.imag;
	temp.imag= lhs.real * rhs.imag + lhs.imag * rhs.real;
	return temp;
}
Complex operator* ( const int &lhs, Complex const & rhs)
{
	Complex temp;
	temp.real=lhs*rhs.real;
	temp.imag=lhs*rhs.imag;
	return temp;
}
istream &operator>>(istream&is, Complex &item)
{
	is >> item.real>>item.imag;
	return is;
}
vector<Complex> recursive_FFT(vector<Complex> const&a)
{
	size_t n = a.size();
	if (n == 1) return a;
	struct  Complex wn, w;
	wn.real = cos(2 * pi / n), wn.imag = sin(2 * pi / n);
	w.real = 1, w.imag = 0;
	vector<Complex>a0, a1;
	a0.reserve(n / 2);
	a1.reserve(n / 2);
	for (int i = 0; i<n / 2; i++)
	{
		a0.push_back(a[i * 2]);
		a1.push_back(a[i * 2 + 1]);
	}
	vector<Complex>y0 = recursive_FFT(a0);
	vector<Complex>y1 = recursive_FFT(a1);
	vector<Complex>y;
	y.resize(n);
	for (int k = 0; k<n / 2; k++)
	{
		y[k] = y0[k] + w*y1[k]; 
		y[k + n / 2] = y0[k] - w*y1[k];
		w = w*wn;
	}
	return y;
}
void main()
{
	cout << "请输入需要计算的多项式最高次数，注意必须是3的幂" << endl;
	int j = 0;
	cin >> j;
	const int n = j;
	vector<Complex>a;
	cout << "请输入多项式各项系数:";
	double  k = 0;
	Complex i;
	while (k != n&&cin >> i)
	{
		a.push_back(i);
		k++;
	}
	vector<Complex>y = recursive_FFT(a);
	for (int t = 0; t<y.size(); t++)
	{
		printf("%g", y[t].real);//求DFT逆需要把结果÷n
		if (fabs(y[t].imag) < 0.001)//虚数小于0.001，那么虚数忽略不计
		{
			cout << endl;
			continue;
		}
		else
		{
			if (y[t].imag<0) printf("%gi", y[t].imag);
			else printf("+%gi", y[t].imag);
		}
		cout << endl;
	}
}

30.2-1证明推论30.4.

30.2-2 计算向量(0,1,2,3)的DFT。

利用上面的FFT程序可得：

30.2-3采用运行时间为O(nlgn)的方案完成练习30.1-1
代码如下：

<span style="color:#000000;">#include <iostream>
#include <stdio.h>
#include <math.h>
#include <vector>
using namespace std;
#define pi 3.1415926535 
struct  Complex
{
	double real;
	double imag;
	int num;
	Complex(int i=0,int j=0){real=i,imag=j;}
};
Complex q(1,0);
Complex operator+ (Complex const & lhs, Complex const & rhs)
{
	Complex temp;
	temp.real =lhs.real+rhs.real;
	temp.imag = lhs.imag+ rhs.imag;
	return temp;
}
Complex operator- (Complex const & lhs, Complex const & rhs)
{
	Complex temp;
	temp.real = lhs.real - rhs.real;
	temp.imag = lhs.imag - rhs.imag;
	return temp;
}
Complex operator* (Complex const & lhs, Complex const & rhs)
{
	Complex temp;
	temp.real = lhs.real * rhs.real - lhs.imag * rhs.imag;
	temp.imag= lhs.real * rhs.imag + lhs.imag * rhs.real;
	return temp;
}
Complex operator* ( const int &lhs, Complex const & rhs)
{
	Complex temp;
	temp.real=lhs*rhs.real;
	temp.imag=lhs*rhs.imag;
	return temp;
}
vector<Complex> operator* (vector<Complex>const&a,vector<Complex>const&b)
{
    vector<Complex> temp;
    for (size_t i=0;i<a.size();i++)
    {
		temp.push_back(a[i]*b[i]);
    }
	return temp;
}
istream &operator>>(istream&is, Complex &item)
{
	is >> item.real>>item.imag;
	return is;
}
vector<Complex> recursive_FFT(vector<Complex> const&a)
{
	size_t n = a.size();//有改动
	if (n == 1) return a;
	struct  Complex wn, w;
	wn.real = cos(2 * pi / n), wn.imag = sin(2 * pi / n);
	w.real = 1, w.imag = 0;
	vector<Complex>a0, a1;
	a0.reserve(n / 2);
	a1.reserve(n / 2);
	for (int i = 0; i<n / 2; i++)
	{
		a0.push_back(a[i * 2]);
		a1.push_back(a[i * 2 + 1]);
	}
	vector<Complex>y0 = recursive_FFT(a0);
	vector<Complex>y1 = recursive_FFT(a1);
	vector<Complex>y;
	y.resize(n);
	for (int k = 0; k<n / 2; k++)
	{
		y[k] = y0[k] + w*y1[k]; 
		y[k + n / 2] = y0[k] - w*y1[k];
		w = w*wn;
	}
	return y;
}
vector<Complex> recursive_FFT_opposite(vector<Complex> const&a)
{
	size_t n = a.size();
	if (n == 1) 
	{
		return a;
	}
	struct  Complex wn, w;
	wn.real = cos(-2 * pi / n), wn.imag = sin(-2 * pi / n);// 注意求DFT的逆是负数的
	w.re