算法导论第三十(30)章多项式与快速傅里叶变换
由于在第三十一章数论算法中遇到几个关于超大数乘法的问题促使我需要学这章，具体请看第三十一章练习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
#include
using namespace std;
struct Polynomial
{
int coef;
//系数
int Index;
//指数
};
int Division_with_remainder(vectora,vector&qx,int x0,const int n)
{
// vectortemp;
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<>x0;
//由于程序并不非常强壮，所以请不要输入字符或者直接输入ctrl+z,只能输入一些整数。输入系数时也是一样的。
cout<<"请按由低次到高次依次输入每项系数，且每项次数均不同(注意系数为0的项不能省略输入，直接补0)："<>x.coef)//输入多项式系数，从低次项到高次项依次输入，输入ctrl+z结束输入。
{
x.Index=i;
a.push_back(x);
i++;
}
const int n=i;
cout<<"您输入的多项式为："<
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
#include
#include
#include
using namespace std;
#define pi 3.1415926535
int m=0;
structComplex
{
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 recursive_FFT(vector const&a)
{
size_t n = a.size();
if (n == 1) return a;
structComplex wn, w;
wn.real = cos(2 * pi / n), wn.imag = sin(2 * pi / n);
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT(a0);
vectory1 = recursive_FFT(a1);
vectory;
y.resize(n);
for (int k = 0;
k> j;
const int n = j;
vectora;
cout << "请输入多项式各项系数:";
doublek = 0;
Complex i;
while (k != n&&cin >> i)
{
a.push_back(i);
k++;
}
vectory = recursive_FFT(a);
for (int t = 0;
t
30.2-1证明推论30.4.
文章图片
30.2-2 计算向量(0,1,2,3)的DFT。
利用上面的FFT程序可得：
文章图片
30.2-3采用运行时间为O(nlgn)的方案完成练习30.1-1
代码如下：
#include
#include
#include
#include
using namespace std;
#define pi 3.1415926535
structComplex
{
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 operator* (vectorconst&a,vectorconst&b)
{
vector temp;
for (size_t i=0;
i>(istream&is, Complex &item)
{
is >> item.real>>item.imag;
return is;
}
vector recursive_FFT(vector const&a)
{
size_t n = a.size();
//有改动
if (n == 1) return a;
structComplex wn, w;
wn.real = cos(2 * pi / n), wn.imag = sin(2 * pi / n);
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT(a0);
vectory1 = recursive_FFT(a1);
vectory;
y.resize(n);
for (int k = 0;
k recursive_FFT_opposite(vector const&a)
{
size_t n = a.size();
if (n == 1)
{
return a;
}
structComplex wn, w;
wn.real = cos(-2 * pi / n), wn.imag = sin(-2 * pi / n);
// 注意求DFT的逆是负数的
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT_opposite(a0);
vectory1 = recursive_FFT_opposite(a1);
vectory;
y.resize(n);
for (int k = 0;
k Polynomial_multiplication(vector const&a,vector const&b)
{
vectorya=recursive_FFT(a);
vectoryb=recursive_FFT(b);
vectoryc=ya*yb;
vectortemp=recursive_FFT_opposite(yc);
return temp;
}
void main()
{
cout << "请输入需要计算的多项式最高次数，注意必须是2的幂" << endl;
int j = 0;
cin >> j;
const int n = j;
vectora,b;
cout << "请输入多项式A各项系数(注意系数为0的项要补0):";
intk1 = 0,k2=0;
Complex i1,i2;
while (k1 != n&&cin >> i1)
{
//x.push_back(i);
a.push_back(i1);
k1++;
}
cout << "请输入多项式B各项系数(注意系数为0的项要补0):";
while (k2 != n&&cin >> i2)
{
b.push_back(i2);
k2++;
}
vectory = Polynomial_multiplication(a,b);
for (int t = 0;
t
结果图：
文章图片
30.2-4 写出伪代码，在O(nlgn)运行时间内计算出DFT^-1.
#include
#include
#include
#include
using namespace std;
#define pi 3.1415926535
structComplex
{
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 recursive_FFT_opposite(vector const&a)
{
size_t n = a.size();
if (n == 1) return a;
structComplex wn, w;
wn.real = cos(-2 * pi / n), wn.imag = sin(-2 * pi / n);
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT_opposite(a0);
vectory1 = recursive_FFT_opposite(a1);
vectory;
y.resize(n);
for (int k = 0;
k> j;
const int n = j;
vectora;
cout << "请输入多项式A各项系数(注意系数为0的项要补0):";
intk = 0;
Complex i;
while (k != n&&cin >> i)
{
a.push_back(i);
k++;
}
vectory = recursive_FFT_opposite(a);
for (int t = 0;
t
30.2-5请把FFT推广到n是3的幂的情形，写出运行时间的递归式并求解。
代码如下：
#include
#include
#include
#include
using namespace std;
#define pi 3.1415926535
structComplex
{
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 recursive_FFT_3(vector const&a)
{
size_t n = a.size();
if (n == 1) return a;
structComplex wn, w,w1,w2;
wn.real = cos(2 * pi / n), wn.imag = sin(2 * pi / n);
//如果求逆需要取负度数
w.real = 1, w.imag = 0;
vectora0,a1,a2;
a0.reserve(n/3);
a1.reserve(n/3);
a2.reserve(n/3);
for (int i = 0;
iy0 = recursive_FFT_3(a0);
vectory1 = recursive_FFT_3(a1);
vectory2 = recursive_FFT_3(a2);
vectory;
y.resize(n);
w1.real=cos(2 * pi / 3),w1.imag=sin(2 * pi /3);
w2.real=cos(4 * pi / 3),w2.imag=sin(4 * pi /3);
for (int k = 0;
k> j;
const int n = j;
vectora,b;
cout << "请输入多项式A各项系数(注意系数为0的项要补0):";
intk1 = 0,k2=0;
Complex i1,i2;
while (k1 != n&&cin >> i1)
{
//x.push_back(i);
a.push_back(i1);
k1++;
}
vectory = recursive_FFT_3(a);
for (int t = 0;
t
文章图片
其递归时间为T(n)=3T(n/3)+O(n)=>T(n)=O(nlgn),所以和基2FFT的复杂度是一样的。
30.2-6略。主要是不太理解”完备的“是什么意思？
30.2-7给定一组值z0,z1...zn-1(可能重复)，说明如何求出仅以z0,z1...zn-1(可能有重复)为零点的一个次数界为n+1的多项式P(x)的系数。你给出的过程运行时间为O(nlgnlgn).(提示：当仅当P(x)在zj处值为0)
由提示可得：P(x)=(x-z0)(x-z1)...(x-zn-1) 将n次单位复数根依次带入多项式，得到向量y=(y1,y2..yn-1)。然后可以用DFT逆求出独享式P(x)的系数。具体看下面代码：
//30.2-7这个版本的时间复杂度是O(nlgnlgn)
#include
#include
#include
#include
using namespace std;
#define pi 3.1415926535
#define nn 8 //输入你想要计算的项数，记住一定要2的幂额
structComplex
{
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 operator* (vectorconst&a,vectorconst&b)
{
Complex y;
vector temp(a.size(),y);
for (size_t i=0;
i operator/(vector&a,double b)
{
vectortemp(a.size(),0);
for (size_t i=0;
i>(istream&is, Complex &item)
{
is >> item.real>>item.imag;
return is;
}
vector recursive_FFT(vector const&a)
{
size_t n = a.size();
//有改动
if (n == 1) return a;
structComplex wn, w;
wn.real = cos(2 * pi / n), wn.imag = sin(2 * pi / n);
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT(a0);
vectory1 = recursive_FFT(a1);
vectory;
y.resize(n);
for (int k = 0;
k recursive_FFT_opposite(vector const&a)
{
size_t n = a.size();
if (n == 1)
{
return a;
}
structComplex wn, w;
wn.real = cos(-2 * pi / n), wn.imag = sin(-2 * pi / n);
// 注意求DFT的逆是负数的
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT_opposite(a0);
vectory1 = recursive_FFT_opposite(a1);
vectory;
y.resize(n);
for (int k = 0;
k Polynomial_multiplication(vector &a,vector &b)
{
Complex y;
int i=0;
int t=a.size();
while (i!=t)
{
a.push_back(y);
b.push_back(y);
i++;
}
vectorya=recursive_FFT(a);
vectoryb=recursive_FFT(b);
vectoryc=ya*yb;
vectortemp=recursive_FFT_opposite(yc);
//应用卷积定理
return temp;
}
vector Polynomial_0(vector const&a)//求(x-z1)(x-z2)..(x-zn)的主函数
{
size_t n=a.size();
if (n==2)//n=2时，可以得到系数向量f(x)=x^2-(z1+z2)x+z1*z2
{
Complex y1(1,0),y2(a[1].real+a[0].real,0),y3(a[1].real*a[0].real,0),y4;
vectoryy(4,y4);
yy[0]=y1;
yy[1]=y2;
yy[2]=y3;
return yy;
}
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0=Polynomial_0(a0);
//分治法分成前一半T(n/2)
vectory1=Polynomial_0(a1);
//分治法分成后一半T(n/2)
vectory;
y.resize(n);
y=Polynomial_multiplication(y0,y1);
//多项式乘法O(nlgn)
y=y/(2*n);
return y;
//T(n)=2T(n/2)+O(nlgn)
}
void main()
{
vectora;
cout << "请输入多项式A各项系数(注意系数为0的项要补0):";
intk1 = 0,k2=0;
Complex i;
while (k1 != nn&&cin >> i)
{
a.push_back(i);
k1++;
}
while (--k1)
{
if (a[k1].real!=0||a[k1].imag!=0)
{
break;
}
k1--;
}
vectory = Polynomial_0(a);
for (int t = 0;
t<=k1+1;
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;
}
}
上面这个版本的时间复杂度是O(nlgn) ，下面的版本是O(n^2)
#include
#include
#include
#include
using namespace std;
#define pi 3.1415926535
int m=0;
structComplex
{
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 recursive_FFT_opposite(vectorconst&a)
{
size_t n = a.size();
if (n == 1)
{
return a;
}
structComplex wn, w;
wn.real = cos(-2 * pi / n), wn.imag = sin(-2 * pi / n);
// 注意求DFT的逆是负数的
w.real = 1, w.imag = 0;
vectora0, a1;
a0.reserve(n / 2);
a1.reserve(n / 2);
for (int i = 0;
iy0 = recursive_FFT_opposite(a0);
vectory1 = recursive_FFT_opposite(a1);
vectory;
y.resize(n);
for (int k = 0;
k recursive_FFT_opposite1(vector const&x,int j)
{//O(n^2)
const int n = j;
vectora;
int k=0;
Complex i,wn,w;
wn.real=cos(2*pi/n),wn.imag=sin(2*pi/n);
w.real=1,w.imag=0;
while (k != n)//O(n)
{//O(n^2)
m=0;
q.real=1,q.imag=0;
for (int t=0;
ty =recursive_FFT_opposite(a);
//O(nlgn)
return y;
}
void main()
{
cout << "请输入需要计算的多项式最高次数，注意必须是3的幂" << endl;
int j = 0;
cin >> j;
const int n = j;
vectora,x;
cout << "请输入多项式A各项系数(注意系数为0的项要补0):";
int k=0;
Complex i,wn,w;
wn.real=cos(2*pi/n),wn.imag=sin(2*pi/n);
w.real=1,w.imag=0;
while (k != n&&cin >> i)
{
x.push_back(i);
//a.push_back(i1);
k++;
}
vectory = recursive_FFT_opposite1(x,j);
for (int t = 0;
t
30.2-8 一个向量a=(a0,a1...an-1)的线性调频变换的向量y=(y0,y1...yn-1),其中yk=∑ajz^(kj),z是任意复数。因此，通过取z=wn,DFT是线性调频变换的一种特殊情形。对任意复数z,请说明如何在O(nlgn)时间内求出线性调频变换的值。 
文章图片
30.3高效?FFT实现 ? ??公用子表达式：多次计算同一个表达式的值，该值就是公用子表达式。 蝴蝶操作：循环中，数据先加上公用子表达式，然后这个数据又减去它，那么这种操作称为蝴蝶操作。 位逆序置换：就是把一个数的二进制位逆序后得到的数。 FFT迭代实现代码如下： //迭代的FFT，比递归的运行更快
#include
#include
#include
using namespace std;
#define pi 3.1415926535
struct Complex
{
double imag;
double real;
Complex(double r=0,double i=0){imag=i,real=r;
}
};
vectora;
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+ (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;
}
int rev(int k)
{
int temp=0,i=0;
static int flag=0;
size_t n=a.size()-1;
if (!flag){while(n>>i)i++;
flag=i;
}//设置哨兵为了只调用一次这个循环
int j=flag;
vectorb(n+1,0);
while (j>0)
{
//b[--j]=(k-((k>>1)<<1));
b[--j]=k%2;
k>>=1;
temp+=b[j]< BIT_REVERSE_COPY(vector&a,vector&A)
{
size_t n=a.size();
for (size_t k=0;
k>=1) temp++;
return temp;
}
vector ITERATIVE_FFT(vector&a,vector&A)
{
BIT_REVERSE_COPY(a,A);
size_t n=a.size();
for (int s=1;
s<=Log(n);
s++)
{
int m=pow(2,s);
for (int k=0;
k<=n-1;
k+=m)
{
Complex w(1,0),wm(cos(2*pi/m),sin(2*pi/m)),t,u;
for (int j=0;
j<=m/2-1;
j++)
{
t=w*A[k+j+m/2];
u=A[k+j];
A[k+j]=u+t;
A[k+j+m/2]=u-t;
w=w*wm;
}
}
}
return A;
}
void main()
{
Complex x,y;
int i=0;
while (cin>>x.real>>x.imag)//输入复系数
{
a.push_back(x);
i++;
}
vectorA(i,y);
BIT_REVERSE_COPY(a,A);
vectoryy = ITERATIVE_FFT(a,A);
for (int t = 0;
t
30.3-1 请说明如何用ITERATIVE-FFT计算出输入向量(0,2,3,-1,4,5,7,9)的DFT。 利用上面FFT迭代法计算出的结果是：
文章图片
30.3-2 请说明如何实现一个FFT算法，注意把位逆序置换放到计算最后而不是在开始。 vector ITERATIVE_FFT(vector&a,vector&A)
{
size_t n=a.size();
for (int s=1;
s<=Log(n);
s++)
{
int m=pow(2,s);
for (int k=0;
k<=n-1;
k+=m)
{
Complex w(1,0),wm(cos(2*pi/m),sin(2*pi/m)),t,u;
for (int j=0;
j<=m/2-1;
j++)
{
t=w*a[rev(k+j+m/2)];
u=a[rev(k+j)];
a[rev(k+j)]=u+t;
a[rev(k+j+m/2)]=u-t;
w=w*wm;
}
}
}
BIT_REVERSE_COPY(a, A);
return A;
}
30.3-3 在每个阶段中，ITERATIVE-FFT计算旋转因子多少次？重写ITERATIVE-FFT,使其在阶段s中计算旋转因子2^(s-1)次。
书上给的迭代FFT，每个阶段需要计算n/2次，重写的代码如下： vector ITERATIVE_FFT(vector&a,vector&A)
{
BIT_REVERSE_COPY(a,A);
size_t n=a.size();
for (int s=1;
s<=Log(n);
s++)
{
int m=pow(2,s);
Complex w(1,0),wm(cos(2*pi/m),sin(2*pi/m)),t,u;
for (int j=0;
j<=m/2-1;
j++)
{
for (int k=0;
k<=n-1;
k+=m)
{
t=w*A[k+j+m/2];
u=A[k+j];
A[k+j]=u+t;
A[k+j+m/2]=u-t;
}
w=w*wm;
}
}
return A;
}
30.3-4 略，没看懂题意。 思考题 30.1(分治乘法)a.说明如何仅用三次乘法，就能求出线性多项式ax+b与cx+d的乘积。(提示：有一个乘法运算时(a+b)(c+d)) (ax+b)(cx+d)=acx^2+((a+b)(c+d)-ac-bd)x+bd b.试写出两种分治算法，求出两个次数界为n的多项式乘积，使其在O(n^lg3)运行时间内。第一个算法把输入多项式系数分成高阶系数一半与低阶系数一半，第二个算法应该根据其系数下标的奇偶性来进行划分。 按高低次数划分情况 //O(n^lg3)的多项式乘法，这个多项式乘法递归函数总体来讲和书上的FFT结构类似
#include
#include
#include
using namespace std;
#define n 8 //n始终是2的幂，不满2的幂的多项式用0为系数的项自动补满使项数达到2的幂
struct Polynomial
{
int coef;
//系数
int index;
//指数
Polynomial(int c=0,int i=0){coef=c,index=i;
};
};
Polynomial operator-(Polynomial&a,Polynomial&b)//运算符重载函数，下面几个函数类似。
{
Polynomial temp;
temp.coef=a.coef-b.coef;
return temp;
}
vector operator+(vector&a,vector&b)
{
vectortemp;
temp.resize(a.size());
for (size_t i=0;
i operator-(vector&a,vector&b)
{
vectortemp;
temp.resize(a.size());
for (size_t i=0;
iPolynomial_multiplication(vector&a,vector&b)//多项式乘法
{//T(n)=3T(n/2)+O(n) =>T(n)=O(n^lg3)
size_t m=a.size();
vectora0,a1,b0,b1;
vectorc;
c.resize(3);
//由于递归式只含3项(f(x)=c[0]+c[1]x^(n/2)+c[2]x^n)，所以c只增加3个元素。
for (size_t s=0;
s>x)
{
a[i].coef=x;
i++;
}
cout<<"请输入项数不超过"<>x)
{
b[j].coef=x;
j++;
}
vector yy=Polynomial_multiplication(a,b);
vectorflag;
flag.resize(yy.size());
int i1=0;
while(i1!=yy.size()) flag[i1++]=y;
for (int k=0;
k
按奇偶次数划分的情况 //O(n^lg3)的多项式乘法(奇偶次数划分)，这个多项式乘法递归函数总体来讲和书上的FFT结构类似。
#include
#include
#include
using namespace std;
#define n 8 //n始终是2的幂，不满2的幂的多项式用0为系数的项自动补满使项数达到2的幂
//int t=-1;
int Log(int nn)
{
int i=0;
while (nn>>=1)i++;
return i;
}
struct Polynomial
{
int coef;
//系数
int index;
//指数
Polynomial(int c=0,int i=0){coef=c,index=i;
};
};
vector operator+(vectorconst &a,vectorconst &b)//运算符重载函数，下面几个函数类似。
{
vectortemp;
temp.resize(a.size());
for ( int i=0;
i operator-(vectorconst &a,vectorconst &b)
{
vectortemp;
temp.resize(a.size());
for (size_t i=0;
iPolynomial_multiplication(vector&a,vector&b)//多项式乘法
{//T(n)=3T(n/2)+O(n) =>T(n)=O(n^lg3)
size_t m=a.size();
vectora0,a1,b0,b1;
vectorc;
c.resize(3);
//由于递归式只含3项(f(x)=c[0]+c[1]x^(n/2)+c[2]x^n)，所以c只增加3个元素。
for (size_t s=0;
s>x)
{
a[i].coef=x;
a[i].index=i;
i++;
}
cout<<"请输入项数不超过"<>x)
{
b[j].coef=x;
b[j].index=j;
j++;
}
vector yy=Polynomial_multiplication(a,b);
cout<
c.证明：请说明如何用O(n^lg3)步计算出两个n位整数的乘积，其中每一步至多常数个1位的值进行操作。 什么叫做“至多常数个1位的值进行操作”？不过下面的算法肯定是O(n^lg3)步计算出两个n位整数的乘积. //O(n^lg3)的多项式乘法转化为两个n位整数乘法，这个多项式乘法递归函数总体来讲和书上的FFT结构类似
#include //更精确的说是(1/2)n^lg3次递归
#include
using namespace std;
#define n 8 //n始终是2的幂，不满2的幂的多项式用0为系数的项自动补满使项数达到2的幂
struct Polynomial
{
int coef;
//系数
int index;
//指数
Polynomial(int c=0,int i=0){coef=c,index=i;
};
};
Polynomial operator-(Polynomial&a,Polynomial&b)//运算符重载函数，下面几个函数类似。
{
Polynomial temp;
temp.coef=a.coef-b.coef;
return temp;
}
vector operator+(vector&a,vector&b)
{
vectortemp;
temp.resize(a.size());
for (size_t i=0;
i operator-(vector&a,vector&b)
{
vectortemp;
temp.resize(a.size());
for (size_t i=0;
iPolynomial_multiplication(vector&a,vector&b)//多项式乘法
{//T(n)=3T(n/2)+O(n) =>T(n)=O(n^lg3)
size_t m=a.size();
vectora0,a1,b0,b1;
vectorc;
c.resize(3);
//由于递归式只含3项(f(x)=c[0]+c[1]x^(n/2)+c[2]x^n)，所以c只增加3个元素。
for (size_t s=0;
s>x)
{
a[i++].coef=x;
}
cout<<"请输入不超过"<>x)
{
b[j++].coef=x;
}
vector yy=Polynomial_multiplication(a,b);
vectorflag;
flag.resize(yy.size());
while(i1!=yy.size()) flag[i1++]=y;
for (int k=0;
k0)//显示整数
{
if(!a[--i].coef&&flag1==0)continue;
else {flag1=1;
cout<0)//显示整数
{
if(!b[--j].coef&&flag1==0)continue;
else {flag1=1;
cout<> j;
const int n = j;
cout<<"打算输入几个有效数据？且不能超过"<>j1;
const int n1=j1;
intk=0;
Complex i2;
vectorb,c,s;
cout << "请输入循环矩阵B(代替向量B)的第一行各项(有效数据"<>i2)
{
b.push_back(i2);
k++;
}
vector >T;
T.resize(n1);
for (size_t i=0;
i>T[0][i];
cout<<"请输入特普利茨矩阵的第一列数据,且数据不能超过"<>T[i][0];
vectory=Toeplitz_vector(T,c,s,b,n1,n);
cout<
d.请给出一个高效算法计算出两个nXn特普利茨矩阵乘积，并分析此算法的运行时间。 利用c中引理和算法规律，计算过程和c中所给代码差不多。只是多计算几个矩阵乘积，但是其本质与c无异。由于求两个矩阵乘积，与c不同的是需要O(n^2)+O(nlgn)次乘法运算。
30-3(多维快速傅里叶变换) 我们可以将式(30.8)定义的一维离散傅里叶变换推广到d维上。这时输入时一个d维的数组A={aj1,j2...jd},维数分别为n1,n2,...nd,其中n1n2...nd=n.定义d维离散傅里叶变换如下：
文章图片
,其中 0≤k1
文章图片
这个元素wn1有n1种数据，wn2有n2种数据。。。。wnd有nd种数据，所以此全排列有n1n2...nd=n种数据，所以向量y一共有n项，其中的每一项yi也是由n个子元素组成。而通项yi中的n个通项子元素
文章图片
是固定的，调整维度顺序仅仅是子元素出现顺序有变化，但是整体的数值是不变的不受维度排列顺序影响，所以通项yi也是固定的，从而得证。 c.证明：如果采用计算快速傅里叶变换计算每个一维的DFT，那么计算一个d维的DFT的总时间是O(nlgn),与d无关。 由a知： (n/n1)O(n1lgn1)+(n/n2)O(n2lgn2)+....(n/nd)O(ndlgnd)=nO(lgn1)+nO(lgn2)+.....+nO(lgnd)=nO(lgn1n2..nd)=O(nlgn)
文章图片
文章图片
文章图片
文章图片
30-5 多项式在多个点的求值 我们已经注意到，运用霍纳法则，就能够在O(n)的时间内，求出次数界为n-1的多项式在单个点的值。同时也发现，运用FFT也能够在O(nlgn)的时间内，求出多项式在所有n个单位复根处的值。现在我们就来说明如何在O(nlg2n)的时间内，求出一个次数界为n的多项式在任意n个点的值。为了做到这一点，我们将不加证明地运用下列结论：当一个多项式除以另一个多项式时，可以在O(nlgn)的时间内计算出其多项式余式。例如，多项式3x3+x2-3x+1除以多项式x2+x+2所得的余式为(3x3+x2-3x+1)mod(x2+x+2)=-7x+5 
文章图片
a)证明 根据带余除法A(x)可以分解为q(x)(x-z)+r(x),所以当x=z时，A(z)=r(z),其中r(z)=A(x)mod(x-z) 所以a式得证！ b)证明 根据Qij=A(x)mod Pij(x) =>Qkk=A(x)mod Pkk(x) =>根据pij(x)=π(x-xk)定义知：Qkk=A(x)mod (x-xk)=A(xk) 的证！Q(0,n-1)=A(x)mod P(0,n-1)(x),其中P(0,n-1)(x)展开后是一个n+1次数界的多项式，而A(x)为次数界为n的多项式，所以A(x) c)证明：根据31.1-7结论 Q(i,j)modP(i,k)=>(A(x)modP(i,j))mod P(i,k)=>(A(x)mod(P(i,k)*P(k,j)))mod P(i,k)=>由于P(i,k) | P(i,j)，所以Q(i,j)modP(i,k)=A(x)mod P(i,k)=Q(i,k),同理另一个等式也可以得证！ d)Q(0,n-1)分解为Q(i,(i+j)/2)与Q((i+j)/2+1,j)，一直递归到Q(i,i),递归式为T(n)=2T(n/2)+O(nlgn) 可以用递归树猜测解并用代入法证明=>T(n)=O(nlg2n),所以所给算法符合题意。 举个例子：/Q(0,0)
/Q(0,1)
\ Q(1,1)
以n=4为例Q(0,3)
/Q(2,2)
\ Q(2,3)
\Q(3,3) 30-6 略。 红字为今天更新的内容，特别地，30.2-7与30.2-8题目更正了错误。 ?? 【算法导论第三十(30)章多项式与快速傅里叶变换】
推荐阅读
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- 画解算法（1.|画解算法：1. 两数之和）
- Guava|Guava RateLimiter与限流算法
- 一个选择排序算法
- SG平滑轨迹算法的原理和实现
- 《算法》-图[有向图]
- LeetCode算法题-11.|LeetCode算法题-11. 盛最多水的容器（Swift）
- 第三十九章|第三十九章 越喜欢的东西
- 遇见第三十四天
- 虚拟DOM-Diff算法详解
- 《数据结构与算法之美》——队列