FFT多项式乘法学习笔记

时间：2015-06-04 09:59:21 阅读：215 评论：0 收藏：0 [点我收藏+]

标签：fft

??其实我不知道我是否真的理解了FFT，但是我会用FFT优化多项式乘法了QAQ。。
（以下大多摘自算导
前置知识
1. 多项式
??在一个代数域F上，关于变量x的多项式定义为形式和形式表示的函数

A (x) = \sum j = 0 n ? 1 a j x j ， 其 中 a 0 \dots a n ? 1 为 多 项 式 各 项 的 系 数

$A(x)=\sum_{j=0}^{n-1}a_jx^j，其中a0…an-1为多项式各项的系数$
2. 多项式的次数界
??若多项式有非零系数的最高次项为

xk $x^k$ ，则称k为该多项式的次数，任何严格大于k的整数都是这个多项式的次数界。
3. 多项式的表示
（1）系数表示法
??对于一个次数界为n的多项式A（x）来说，其系数表示法可以看做是一个列向量

a=（a0，a1，…，an?1） $a=（a_0，a_1，…，a_{n-1}）$ 。系数表示法对于某些多项式的计算很方便，如对多项式A（x）在给定点

x0 $x_0$ 的求值运算就是计算A（x0）的值。如果使用霍纳法则，则求值运算的运行时间为

O（n） $O（n）$ ，

A （ x 0 ） = a 0 + x 0 (a 1 + x 0 (a 2 + . . . + x 0 (a n ? 2 + x 0 (a n ? 1)) . . .))

$A（x_0）=a_0+x_0(a_1+x_0(a_2+...+x_0(a_{n-2}+x_0(a_{n-1}))...))$
??另外，加法运算的时间复杂度是

O（n） $O（n）$ ，暴力进行乘法运算的时间复杂度是

O（n2） $O（n^2）$ 。
（2）点值表示法
?? 有n个点

（x0，A（x0），），（x1，A（x1）），..，（xn?1，A（n?1）） ${（x_0，A（x_0），），（x_1，A（x_1）），..，（x_{n-1}，A（n-1））}$ ，当所有

xk $x_k$ 各不相同时，这n个点可以唯一表示一个次数界为n的多项式，但是一个次数界为n的多项式可以有多个点值表示。通俗一点说，已知一个多项式函数的n个函数值可以唯一确定这个函数，而知道这个函数可以知道不止n个函数值。
??已知点值表示求系数表示称为插值，用拉格朗日插值法可以做到

O（n2） $O（n^2）$ 。
??若次数界为n的多项式A（x），B（x）的n个点的

xk $x_k$ 是对应相同的，点值表示法的加法操作时间复杂度是O（n），只要把对应A/B（

xk $x_k$ ）相加即可，若A（x），B（x）都已知

2n $2n$ 个点，那把A/B（

xk $x_k$ ）相乘同样可以在

O（n） $O（n）$ 时间内完成乘法运算。。

进入正题。。
1. 单位复根
??n次单位复根是满足

wn=1 $w^n=1$ 的复数

w $w$ ，有n个，他们均匀的分布在以复
平面的原点为圆心的单位圆上，为

e2πki/n，k=0,1，…，n?1 $e^{2πki/n}，k=0,1，…，n-1$ ，复数幂定义为

eui=cos(u)+sin(u)i $e^{ui}=cos(u)+sin(u)i$ 。

wn=e2πi/n $w_n=e^{2πi/n}$ 称为主n次单位根，所有的其他n次单位根都是

wn $w_n$ 的幂。
??因为有

wnn=w0n=1 $w_n^n=w_n^0=1$ ，所以有

wjnwkn=w(j+k)modnn $w_n^jw_n^k=w_n^{(j+k)mod n}$ 。（性质1

重要性质
相消引理
?? 对任何整数

n≥0,k≥0,d>0 $n\ge0,k\ge0,d>0$ ，有

wdkdn=wkn $w_{dn}^{dk}=w_n^k$ ，利用定义不难证明
?? 推论：对任意偶数

n>0 $n>0$ ，有

wn/2n=w2=1 $w_n^{n/2}=w_2=1$ 。
折半引理
?? 如果

n>0 $n>0$ 为偶数，n个n次单位复根的平方等于

n/2 $n/2$ 个

n/2 $n/2$ 次单位复根。利用相消引理和性质1可证。
求和引理
??对于任意整数

n≥1 $n\ge1$ 和不能被n整除的非零整数k，有

∑n?1j=0wkjn=0 $\sum_{j=0}^{n-1}w_n^{kj}=0$ 。
证明：
?? 原式

=(wkn)n?1wkn?1=(wnn)k?1wkn?1=1?1wkn?1=0 $=\frac{(w_n^k)^n-1}{w_n^k-1}=\frac{(w_n^n)^k-1}{w_n^k-1}=\frac{1-1}{w_n^k-1}=0$ ，只要k不整除n就可以保证分母不为0。

2. DFT
??A（x）是一个次数界为n的多项式，不失一般性地假定n是2的幂，因为
次数界总是可以增大的。有列向量

y=（y0，y1，…，yn?1） $y=（y_0，y_1，…，y_{n-1}）$ ，其中

yk=A（wkn） $y_k=A（ w_n^k）$ ，则称y是系数向量a的离散傅里叶变换，也写作

y=DFTn(a) $y=DFT_n (\bf{a})$ 。

3. 快速傅里叶变换（FFT）
??FFT是一种可以在

O（nlogn） $O（nlogn）$ 时间内计算出

DFTn(a) $DFT_n (\bf{a})$ 的算法，主要思想
是分治。
??定义

A0（x）=a0+a2x+a4x2+…+an?2xn/2?1 $A_0（x）=a_0+a_2x+a_4x^2+…+a_{n-2}x^{n/2-1}$ ，包含了A（x）
所有偶数下标的系数，

A1（x）=a1+a3x+a5x2+…+an?1xn/2?1 $A_1（x）= a_1+a_3x+a_5x^2+…+a_{n-1}x^{n/2-1}$ ,包含了A（x）所有奇数下标的系数。易得

A（x）=A0（x2）+xA1（x2） $A（x）=A_0（x^2）+xA_1（x^2）$ ，所以我们可以先求

A0 $A_0$ 和

A1 $A_1$ 的DFT，然后再组合起来。组合的时候有

A (w k + n / 2 n) = A 0 (w 2 k + n n) + w k + n / 2 n A 1 (w 2 k + n n) = A 0 (w k n) + w k n w 2 A 1 (w k n) A 0 (w k n) ? w k n A 1 (w k n)

$A(w_n^{k+n/2})=A_0(w_n^{2k+n})+w_n^{k+n/2}A_1(w_n^{2k+n})=A_0(w_n^k)+w_n^kw_2A_1(w_n^k)A_0(w_n^k)-w_n^kA_1(w_n^k)$
所以用

n2 $\frac{n}{2}$ 个

yk $y_k$ 就可以推出全部。根据折半引理，问题的规模缩小一半，每次组合的时间复杂度是

O（n） $O（n）$ ，所以总时间复杂度是

O（nlogn） $O（n log n）$ 。
?? 求出

DFTn(a) $DFT_n (\bf{a})$ 后，要对单位复根进行插值，将点值表示转化为系数表示，有

y=Vna $\bf{y}=V_n\bf{a}$ ,其中

V n = ? ? ? ? ? ? ? ? ? ? ? ? 1111 . . . 1 1 w n w 2 n w 3 n . . . w n ? 1 n 1 w 2 n w 4 n w 6 n . . . w 2 n 1 w 3 n w 6 n w 9 n . . . w 3 n . . . . . . . . . . . . . . . . . . 1 w n ? 1 n w 2 (n ? 1) n w 3 (n ? 1) n . . . w (n ? 1) (n ? 1) n ? ? ? ? ? ? ? ? ? ? ? ?

$V_n= \begin{equation} \left[ \begin{array}{ccc} 1&1&1&1&...&1\ 1&w_n&w_n^2&w_n^3&...&w_n^{n-1}\ 1&w_n^2&w_n^4&w_n^6&...&w_n^{2(n-1)}\ 1&w_n^3&w_n^6&w_n^9&...&w_n^{3(n-1)}\ ...&...&...&...&...&...\1&w_n^{n-1}&w_n^2&w_n^3&...&w_n^{(n-1)(n-1)} \end{array} \right] \end{equation}$
是范德蒙特矩阵。

a=Vn?1?y $\bf{a}=V_{n-1}*\bf{y}$ ，考虑

Vn?1 $V_{n-1}$ 在

（j，k） $（j，k）$ 处的数

Vj,k $V_{j,k}$ ，根据

V?1nVn=In $V_n^{-1}V_n=I_n$ （n阶单位矩阵）有

\sum t = 0 n ? 1 V j, t w k t n = {1 (j = k) 0 (j! = k)

$\sum_{t=0}^{n-1}V_{j,t}w_n^{kt}= \begin{equation} \left\{ \begin{array}{cc} 1(j=k)\ 0(j!=k) \end{array} \right . \end{equation}$
定理：

Vj,k=w?kjn/n $V_{j,k}=w_n^{-kj}/n$
证明：

∑n?1t=0Vj,twktn=1n∑n?1t=0wt(k?j)n $\sum_{t=0}^{n-1}V_{j,t}w_n^{kt}=\frac{1}{n}\sum_{t=0}^{n-1}w_n^{t(k-j)}$ ,当

j=k $j=k$ 时，上式=1，否则由求和引理可得上式为0，得证。
??那么有

aj=1n∑n?1k=0ykw?kjn $a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_kw_n^{-kj}$ ，因此在求逆DFT的时候可以类似于求DFT，只要把y和a角色互换，然后让

wn=w?1n $w_n=w_n^{-1}$ ，再做FFT即可。
??写的时候发现非递归要比递归快很多。。

递归：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<memory.h>
#define N 400010
using namespace std;
const double pi=acos(-1);
struct complex{
    double x,i;
    complex(){}
    complex(double x,double i):x(x),i(i){}
    complex operator+(complex a) {return complex(x+a.x,i+a.i);}
    complex operator-(complex a) {return complex(x-a.x,i-a.i);}
    complex operator*(complex a) {return complex(x*a.x-i*a.i,x*a.i+i*a.x);}
}a[N],b[N];
int n,m,i,nn;
void fft(complex *a,int n,int t)
{
    if (n==1) return;
    complex a0[n>>1],a1[n>>1];
    for (int i=0;i<n;i+=2) a0[i>>1]=a[i],a1[i>>1]=a[i+1];
    fft(a0,n>>1,t);fft(a1,n>>1,t);
    complex wn(cos(2*pi/n),t*sin(2*pi/n)),w(1,0);
    for (int i=0;i<(n>>1);i++,w=w*wn) a[i]=a0[i]+w*a1[i],a[i+(n>>1)]=a0[i]-w*a1[i];
}
int main()
{
    scanf("%d%d",&n,&m);
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
    for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
    nn=1;while (nn<=n+m) nn<<=1;
    fft(a,nn,1);fft(b,nn,1);
    for (i=0;i<=nn;i++) a[i]=a[i]*b[i];
    fft(a,nn,-1);
    for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/nn+0.5));
}

非递归：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<memory.h>
#define N 400010
using namespace std;
const double pi=acos(-1);
struct complex{
    double x,i;
    complex(){}
    complex(double x,double i):x(x),i(i){}
    complex operator+(complex a) {return complex(x+a.x,i+a.i);}
    complex operator-(complex a) {return complex(x-a.x,i-a.i);}
    complex operator*(complex a) {return complex(x*a.x-i*a.i,x*a.i+i*a.x);}
}a[N],b[N];
int n,m,i,nn,len,rev[N];
void fft(complex *a,int n,int t)
{
    for (int i=0;i<n;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int j=1;j<n;j<<=1)
    {
        complex wn(cos(2*pi/(j<<1)),t*sin(2*pi/(j<<1)));
        for (int i=0;i<n;i+=(j<<1))
        {
            complex w(1,0),t0,t1;
            for (int k=0;k<j;k++,w=w*wn) t0=a[i+k],t1=w*a[i+j+k],a[i+k]=t0+t1,a[i+j+k]=t0-t1;
        }
    }
}
int main()
{
    freopen("FFT.in","r",stdin);
    scanf("%d%d",&n,&m);
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for (i=0;i<=n;i++) scanf("%lf",&a[i].x);
    for (i=0;i<=m;i++) scanf("%lf",&b[i].x);
    nn=1;len=0;while (nn<=n+m) nn<<=1,len++;
    rev[0]=0;
    for (i=1;i<nn;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
    fft(a,nn,1);fft(b,nn,1);
    for (i=0;i<=nn;i++) a[i]=a[i]*b[i];
    fft(a,nn,-1);
    for (i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/nn+0.5));
}

??非递归这里有一个翻转的函数，意在让所有数按合并时候的二叉树的叶子节点的顺序排列，不难发现翻转过来就是它的新位置。。

FFT多项式乘法学习笔记

标签：fft

原文地址：http://blog.csdn.net/tag_king/article/details/46351821

踩

(0)

评论一句话评论（0）

分享档案

更多>

2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)

周排行