# 多项式乘法优化 学习笔记

①如何用一次DFT加一次IDFT求出两个实序列A和B的卷积？

${C}^{2}\left[k\right]=\sum _{j=0}^{k}C\left[j\right]C\left[k-j\right]$

$=\sum _{j=0}^{k}\left(A\left[j\right]+i\ast B\left[j\right]\right)\left(A\left[k-j\right]+i\ast B\left[k-j\right]\right)$

$=\sum _{j=0}^{k}\left(A\left[j\right]A\left[k-j\right]-B\left[j\right]B\left[k-j\right]\right)+i\ast \sum _{j=0}^{k}\left(A\left[j\right]B\left[k-j\right]+B\left[j\right]A\left[k-j\right]\right)$

②如何用一次DFT同时求出两个实序列在单位复数根处的点值？

$P\left(x\right)=A\left(x\right)+i\ast B\left(x\right)$

$Q\left(x\right)=A\left(x\right)-i\ast B\left(x\right)$

$P\left(x\right),Q\left(x\right)$$P(x),Q(x)$在单位复数根处的点值表达分别为${F}_{P},{F}_{Q}$$F_P,F_Q$，则可以证明${F}_{P}\left(k\right)$$F_P(k)$${F}_{Q}\left(N-k\right)$$F_Q(N-k)$互为共轭复数。因此只需要对$P\left(x\right)$$P(x)$进行DFT即可。然后会有：

$DF{T}_{A}=\frac{{F}_{P}+{F}_{Q}}{2}$

$DF{T}_{B}=\frac{{F}_{P}-{F}_{Q}}{2i}$

（myy论文里第二条式子写的是乘以$i$$i$，不过我觉得是除以$i$$i$才对）

③关于单位复数根${\omega }_{N}^{k}$$\omega_N^k$

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
#include<vector>
using namespace std;

const int maxn=4000000;
const double pi=acos(-1.0);

struct Complex
{
double X,Y;
Complex (double a=0.0,double b=0.0) : X(a),Y(b) {}
} ;

Complex operator+(Complex a,Complex b){return Complex(a.X+b.X,a.Y+b.Y);}
Complex operator-(Complex a,Complex b){return Complex(a.X-b.X,a.Y-b.Y);}
Complex operator*(Complex a,Complex b){return Complex(a.X*b.X-a.Y*b.Y,a.X*b.Y+a.Y*b.X);}

Complex A[maxn];
Complex B[maxn];

vector <Complex> w[maxn];
int Rev[maxn];
int N,Lg;

int F[maxn];
int G[maxn];
int n,m;

void DFT(Complex *a,double f)
{
for (int i=0; i<N; i++)
if (i<Rev[i]) swap(a[i],a[ Rev[i] ]);

for (int len=2; len<=N; len<<=1)
{
int mid=(len>>1);
for (Complex *p=a; p!=a+N; p+=len)
for (int i=0; i<mid; i++)
{
Complex temp=w[mid][i];
if (f==-1.0) temp.Y=-temp.Y;
temp=temp*p[mid+i];
p[mid+i]=p[i]-temp;
p[i]=p[i]+temp;
}
}
}

void FFT()
{
N=1,Lg=0;
while (N<n+m+4) N<<=1,Lg++;
for (int i=0; i<N; i++)
for (int j=0; j<Lg; j++)
if (i&(1<<j)) Rev[i]|=(1<<(Lg-j-1));

int len=1;
while ((len<<1)<=N)
{
double ang=pi/len;
for (int i=0; i<len; i++)
w[len].push_back( Complex( cos(ang*(double)i) , sin(ang*(double)i) ) );
len<<=1;
}

for (int i=0; i<N; i++) A[i]=Complex((double)F[i],(double)G[i]);
DFT(A,1.0);
for (int i=0; i<N; i++) B[i]=A[(N-i)%N],B[i].Y=-B[i].Y;
for (int i=0; i<N; i++)
{
Complex a=A[i]+B[i];
a.X/=2.0;
a.Y/=2.0;
B[i]=A[i]-B[i];
B[i].X/=2.0;
B[i].Y/=2.0;
swap(B[i].X,B[i].Y);
B[i].Y=-B[i].Y;
A[i]=a;
}
for (int i=0; i<N; i++) A[i]=A[i]*B[i];
DFT(A,-1.0);
for (int i=0; i<N; i++) A[i].X/=((double)N);
for (int i=0; i<N; i++) F[i]=(int)floor( A[i].X+0.5 );
}

int main()
{
freopen("3803.in","r",stdin);
freopen("3803.out","w",stdout);

scanf("%d%d",&n,&m);
for (int i=0; i<=n; i++) scanf("%d",&F[i]);
for (int i=0; i<=m; i++) scanf("%d",&G[i]);

FFT();
for (int i=0; i<=n+m; i++) printf("%d ",F[i]);
printf("\n");

return 0;
}

©️2019 CSDN 皮肤主题: 技术黑板 设计师: CSDN官方博客