zcmimi's blog

FFT

fst fst tle

DFT: 离散傅里叶变换

IDFT: 离散傅里叶逆变换

FFT: 快速傅里叶变换

FNTT/NTT: 快速傅里叶变换的优化版

FWT: 快速沃尔什变换,利用类似FFT的东西解决一类卷积问题

MTT: 毛爷爷的FFT,非常nb/任意模数

FMT: 快速莫比乌斯变化

(摘自https://www.cnblogs.com/zwfymqz/p/8244902.html)

为什么要用到FFT呢?

以高精度乘法举个例子:

你现在要计算a\times b,a,b>10^{1000000}

len_a=n,len_b=m

如果用普通的高精度乘法\mathcal{O(nm)},直接tle

但是FFT可以做到\mathcal{O((n+m) \log (n+m))}

原理: 先把多项式转化为用点值表示\mathcal{O(n \log n)}

然后再用点值相乘来计算\mathcal{O(n)}

然后再通过点值还原多项式

下面是前置知识:

多项式

系数表示法

形式: f(x)=\sum_{i=0}^{n-1} a_ix^i

比如123456也可以看做一个多项式:

f(10)=6+5\cdot10+4\cdot10^2+3\cdot10^3+2\cdot 10^4+1\cdot 10^5

点值表示法

n个点x_1,x_2, \dots ,x_n

代入f(x)会得到n个不同的y_i

f(x)(x_1,y_1),(x_2,y_2),\dots (x_n,y_n)唯一确定

定理:

一个n-1次多项式在n个不同点的取值唯一确定了该多项式。

证明: (反证法)

假设命题不成立,那么存在两个不同的n-1次多项式A(x),B(x)满足: \forall i\in[1,n], A(x_i)=B(x_i)

C(x)=A(x)-B(x),那么C(x)也是一个n-1次多项式,且\forall i\in[1,n], C(x_i)=0

这样的话C(x)n个根,与代数基本定理(n次多项式有且仅有n个根)不符,所以C(x)不是n-1次多项式

所以原命题成立

(摘自https://zhuanlan.zhihu.com/p/31584464)

复数

复数由实部(real)和虚部(image)组成,形如a+bi

i是虚树单位,i^2=-1,i=\sqrt{-1}

复数运算法则:

(a+bi)+(c+di)=(a+c)+(b+d)i

(a+bi)-(c+di)=(a-c)+(b-d)i

(a+bi)\cdot(c+di)=(ac-bd)+(ad+bc)i

复数z=a+bi在复平面中对应坐标为(a,b)

复平面中x轴又称实轴,y轴又称虚轴

共轭复数

复数z=a+bi的共轭复数z'=a-bi,z'称为z的复共轭

性质:

z\cdot z' = a^2 + b^2

|z|=|z'|(模长相等,到原点的距离相等)

单位根

\omega^n=1\omegan次单位根

在复平面上,以原点为圆心,1为半径作圆,所得的圆叫单位圆。

以圆心为起点,圆的n等分点为终点,做n个向量。

其中幅角为正且最小的向量称为n次单位向量,记为\omega_n^1

根据复数乘法的运算法则,其余n−1个复数为\omega_n^2,\omega_n^3,\dots,\omega_n^n

\omega_n^0=\omega_n^n=1,在复平面上对应的点为(1,0)

–1.5–1.5–1.5–1–1–1–0.5–0.5–0.50.50.50.51111.51.51.52222.52.52.5–1–1–1–0.5–0.5–0.50.50.50.5111000θθθOOOωωω888(cos θ, i sin θ)(cos θ, i sin θ)(cos θ, sin θ)CCCBBBDDDEEEFFFGGG

向量\overrightarrow{O\omega_8}表示的复数就是8次单位根,\omega_8^1

拓展:

欧拉公式: e^{i\theta}=\cos\theta + i \sin \theta(在复平面内坐标为(\cos \theta, \sin \theta))

\omega_n^k=e^{i 2\pi \frac kn} (可以理解为从\omega_n^0(1,0)开始走长度为2\pi \frac kn的圆弧后到达的位置)

欧拉公式: e^{i\pi}=-1(复平面中的(-1,0))

关于证明和理解,可以参考:

数学常数与欧拉公式

https://www.bilibili.com/video/av63666593

https://www.bilibili.com/video/av79134103

性质:

\omega_n^k= \cos 2\pi\frac kn+ i \sin 2\pi \frac kn

\omega_{2n}^{2k}=\omega_n^k

\omega_{2n}^{2k}=\cos 2\pi\frac {2k}{2n}+ i \sin 2\pi \frac {2k}{2n}=\cos 2\pi\frac kn+ i \sin 2\pi \frac kn=\omega_n^k

\omega_n^0=\omega_n^n=1

\omega_n^{k+\frac n2}=-\omega_n^k

\omega _{n}^{\frac{n}{2}}=\cos\frac{n}{2}\cdot\frac{2\pi}{n}+i\sin\frac{n}{2}\cdot\frac{2\pi}{n}=-1 (也就是上面说的e^{i\pi})

正题开始:

DFT

将多项式从系数表示法转换为点值表示法

如果暴力代入n个点,复杂度还是\mathcal{O(n^2)}

我们可以利用单位根的性质

假设设n为偶数

对于A(x)=a_0+a_1x+a_2x^2+a_3x^3+\dots+a_{n-1}x^{n-1}

A(x)=(a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+...+a_{n-1}x^{n-1})

A_1(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^\frac{n-2}2,A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^\frac{n-2}2

那么A(x)=A_1(x^2)+xA_2(x^2)

0 \le k \le \frac{n}{2}-1,k\in Z

A(\omega_n^k)=A_1(\omega_\frac n2 ^k)+\omega_n^k\cdot A_2(\omega_\frac n2 ^k)

对于\frac n2 \le k+\frac n2 \le n-1

A(\omega_n^{k+\frac n2})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac n2}\cdot A_2(\omega_n^{2k+n})

\because \omega_n^{k+\frac n2}=-\omega_n^k,\omega_n^{2k+n}=\omega_\frac n2^k

\therefore A(\omega_n^{k+\frac n2})=A_1(\omega_\frac n2^k)-\omega_n^k\cdot A_2(\omega_\frac n2^k)

如果已经知道A_1,A_2分别在\omega_\frac n2 ^{0,1,2,...,\frac n2 -1}的取值,那么就可以在\mathcal{O(n)}的时间内计算出A(x),而A_1(x),A_2(x)的规模都是A(x)的一半

复杂度\mathcal{T(n)=2T(\frac n2)+O(n) = O(n \log n)}

IDFT

使用快速傅里叶变换将点值表示的多项式转化为系数表示,这个过程叫做离散傅里叶逆变换

即由n维点值向量(A(x_0),A(x_1),\dots,A(x_{n-1}))推出n维系数向量(a_0,a_1,\dots,a_{n-1})

y_0,y_1,y_2,...,y_{n-1}是多项式A(x)变换后的点值

(d_0,d_1,\dots,d_{n-1})(a_0,a_1,\dots,a_{n-1})离散傅里叶变换后的结果

构造F(x)=d_0+d_1x+d_2x^2+\dots+d_{n-1}x^{n-1}

设向量(c_0,c_1,\dots,c_{n-1})

c_kF(x)x=\omega_n^{-k}时的点值表示

c_k=\sum_{i=0}^{n-1}d_i\cdot(\omega_n^{-k})^i \\ =\sum_{i=0}^{n-1}[\sum_{j=0}^{n-1}a_j(\omega_n^i)^j]\cdot(\omega_n^{-k})^i \\ =\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^i)^j\cdot(\omega_n^{-k})^i \\ =\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}

观察\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}

j=k\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}=\sum_{i=0}^{n-1}1=n

j\not = k

\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}=

\sum_{i=0}^{\frac{n-1}2}(\omega_n^i)^{j-k}+(\omega_n^{i+\frac{n-1}2})^{j-k}=0

也就是n次单位根会互相抵消

\therefore \sum_{i=0}^{n-1}(\omega_n^i)^{j-k}=[j=k]\cdot n

\therefore c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^i)^{j-k}=\sum_{j=0}^{n-1}a_j [j=k]\cdot n=na_k

\therefore a_k=\frac{c_k}n

那么我们要将点值表示的多项式转化为系数表示:

  1. (\omega_n^0,\omega_n^{-1},\omega_n^{-2},\dots,w_n^{-(n-1)})作为差值节点,将点值再做一次DFT

    得到(c_0,c_1,c_2,\dots,c_{n-1})

  2. (\frac{c_0}n,\frac{c_1}n,\frac{c_2}n,\dots,\frac{c_{n-1}}n)就是原来的(a_0,a_1,a_2,\dots,a_{n-1})

到这里我们就已经大致懂得FFT的原理了,那要怎么实现呢?

实现

LG 3803 【模板】多项式乘法(FFT)

普通的递归写法

FstFstTLE

不断对下标进行奇偶分类,分成两个子序列之后再不断合并

需要辅助数组,容易TLE,MLE

迭代实现

可以发现在不断进行奇偶分类之后,原数在数组中的位置下标最终变成了下标的二进制反转

我们可以直接通过枚举来代替递归

从最底层开始向上合并,那么怎么合并呢?

蝴蝶操作

设该层一共有n项需要处理

A_1(\omega_{\frac n2}^k)A_2(\omega_{\frac n2}^k)分别存放在a_ka_{k+\frac n2}

A_1(\omega_n^k)A_2(\omega_n^k)要存放在buf(k)buf(k+\frac n2)

我们只需设中间变量t=\omega_n^k\cdot a(k+\frac n2)

buf(k+\frac n2)= a(k) - t

buf(k) = a(k)+t

我们可以发现我们不需要辅助存放的数组buf了,直接在原序列操作即可

代码:

#include<bits/stdc++.h>
#define db double
#define For(i,x,y) for(int i=x;i<=y;++i)
const int N=4000011;
const db Pi=acos(-1.0);
int n,m,limit=1,l=0,r[N];
struct cp{db r,i;cp(db R=0,db I=0){r=R,i=I;}}a[N],b[N];
cp operator + (cp a,cp b){return cp(a.r+b.r,a.i+b.i);}
cp operator - (cp a,cp b){return cp(a.r-b.r,a.i-b.i);}
cp operator * (cp a,cp b){return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
void fft(cp *A,int typ){
    For(i,0,limit-1)
        if(i<r[i])std::swap(A[i],A[r[i]]);
    for(int len=1;len<limit;len<<=1){//已处理好的长度块,也就是n/2
        cp Wn(cos(Pi/len),typ*sin(Pi/len));//单位根
        for(int i=0;i<limit;i+=(len<<1)){
            cp w(1,0);//w_n^0
            For(k,0,len-1){//蝴蝶操作
                cp t=w*A[i+k+len];
                A[i+k+len]=A[i+k]-t;
                A[i+k]=A[i+k]+t;
                w=w*Wn;
            }
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    For(i,0,n)scanf("%lf",&a[i].r);
    For(i,0,m)scanf("%lf",&b[i].r);

    while(limit<=n+m)limit<<=1,++l;
    For(i,0,limit-1)//预处理出不断奇偶分类后的最终位置
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

    fft(a,1);fft(b,1);//DFT
    For(i,0,limit)a[i]=a[i]*b[i];//点值相乘
    fft(a,-1);//IDFT
    For(i,0,n+m)printf("%d ",int(a[i].r/limit+0.5));
}

本文部分内容摘自:(参考资料)

https://zhuanlan.zhihu.com/p/31584464

https://www.cnblogs.com/zwfymqz/p/8244902.html

快速傅里叶变换 FFT
comment评论
Search
search