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

### DFT

CLRS上提供了一种很好的方式，也是常用方法之一

此时$\omega_n=e^{\frac{2i\pi}{n}}$

$1.\omega_n^n=1$

$2.\omega_n^k$互不相等

$3.\omega_n^2=\omega_{\frac{n}{2}},\omega_n^{\frac{n}{2}+k}=-\omega_n^k$

4.\sum_{k=0}^{n-1}(\omega_n^{j-i})^k=\begin{eqnarray*}\left\{\begin{aligned}0,~~~&i\neqj\\n,~~~&i=j\end{aligned}\right.\end{eqnarray*}

$A(\omega_n^k)=$

$T(n)=2T(\frac{n}{2})+O(n)=O(nlogn)$

void rev(int n,cpx*t){
for(int i=0,j=0;i!=n;i++){
if(i>j)swap(t[i],t[j]);
for(int l=n/2;(j^=l)<l;l>>=1);
}
}
void dft(int n,cpx*x,cpx*w){
rev(n,x);
for(int i=2;i<=n;i<<=1){
int m=i/2;
for(int j=0;j<n;j+=i){
for(int k=0;k!=m;k++){
cpx t=x[j+m+k]*w[n/i*k];
x[j+m+k]=x[j+k]-t;
x[j+k]=x[j+k]+t;
}
}
}
}

struct cpx {
double r,i;
cpx(double real=0.0,double image=0.0) {
r=real;i=image;
}
cpx operator +(const cpx w) {
return cpx(r+w.r,i+w.i);
}
cpx operator -(const cpx w) {
return cpx(r-w.r,i-w.i);
}
cpx operator *(const cpx w) {
return cpx(r*w.r-i*w.i,r*w.i+i*w.r);
}
};

### 多项式乘法

#include<map>
#include<stack>
#include<queue>
#include<cstdio>
#include<string>
#include<vector>
#include<cstring>
#include<complex>
#include<iostream>
#include<assert.h>
#include<algorithm>
using namespace std;
#define inf 1001001001
#define infll 1001001001001001001LL
#define ll long long
#define dbg(vari) cerr<<#vari<<" = "<<(vari)<<endl
#define gmax(a,b) (a)=max((a),(b))
#define gmin(a,b) (a)=min((a),(b))
#define Ri register int
#define gc getchar()
#define il inline
bool f=true;Ri x=0;char ch;while(!isdigit(ch=gc))if(ch=='-')f=false;while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=gc;}return f?x:-x;
}
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout);
#define N 888888
struct cpx {
double r,i;
cpx(double real=0.0,double image=0.0) {
r=real;i=image;
}
cpx operator +(const cpx w) {
return cpx(r+w.r,i+w.i);
}
cpx operator -(const cpx w) {
return cpx(r-w.r,i-w.i);
}
cpx operator *(const cpx w) {
return cpx(r*w.r-i*w.i,r*w.i+i*w.r);
}
};
cpx a[N],b[N];
cpx epsilon[N],arti_epsilon[N];
void init_epsilon(int n){
double pi=acos(-1);
for(int i=0;i!=n;i++){
epsilon[i]=cpx(cos(2.0*pi*i/n),sin(2.0*pi*i/n));
//arti_epsilon[i]=conj(epsilon[i]);
arti_epsilon[i]=cpx(cos(2.0*pi*i/n),-sin(2.0*pi*i/n));
}
}
void rev(int n,cpx*t){
for(int i=0,j=0;i!=n;i++){
if(i>j)swap(t[i],t[j]);
for(int l=n/2;(j^=l)<l;l>>=1);
}
}
void dft(int n,cpx*x,cpx*w){
rev(n,x);
for(int i=2;i<=n;i<<=1){
int m=i/2;
for(int j=0;j<n;j+=i){
for(int k=0;k!=m;k++){
cpx t=x[j+m+k]*w[n/i*k];
x[j+m+k]=x[j+k]-t;
x[j+k]=x[j+k]+t;
}
}
}
}
cpx c[N];
void mul(cpx *a,cpx *b){
int A,B;
A=gi;B=gi;++A,++B;
for(int i=0;i<A;i++)a[i].r=gi;
for(int i=0;i<B;i++)b[i].r=gi;
int t=max(A,B);
int n=1;
for(;n<t*2;n<<=1);
init_epsilon(n);
dft(n,a,epsilon);
dft(n,b,epsilon);
for(int i=0;i<n;i++)c[i]=a[i]*b[i];
int nn=n;cout<<n;
dft(n,c,arti_epsilon);
for(int i=0;i<=A+B-1;i++)printf("%d ",(int)(c[i].r/nn+0.5));
}
int main(){
mul(a,b);
return 0;
}

### FNT

posted @ 2016-12-21 20:29 zhouyis 阅读(...) 评论(...) 编辑 收藏