FFT感性瞎扯

为了自己以后再用FFT时不再一脸懵X，本蒟蒻决定感性理解一下FFT。

FFT可以干啥？

把两个多项式乘在一起。

具体地说，对于两个多项式

f : f (x) = Σ_{i = 0}^{n} f_{i} x^{i}

$f:f(x)=\Sigma_{i=0}^n f_i x^i$ 和

g : g (x) = Σ_{i = 0}^{m} g_{i} x^{i}

$g:g(x)=\Sigma_{i=0}^m g_i x^i$ ，得到一个多项式

h : h (x) = Σ_{x = 0}^{n + m} Σ_{i = 0}^{x} f (i) g (x - i)

$h:h(x)=\Sigma_{x=0}^{n+m}\Sigma_{i=0}^x f(i)g(x-i)$ 。

显然，如果暴力的话，是

O (n^{2})

$O(n^2)$ 的。

但是，FFT可以做到 $O(nlog_2n)$ ！

我们一点一点把它扯个不明不白。

Part1.复数

只需要记住一点结论：

复数相乘时，模长相乘，辐角相加。

比如说，对于两个极角表示法的复数

u (u_{l}, u_{θ})

$u(u_l,u_{\theta})$ 与

v (v_{l}, v_{θ})

$v(v_l,v_{\theta})$ ，有

w = u v = (u_{l} v_{l}, u_{θ} + v_{θ})

$w=uv=(u_l v_l,u_{\theta}+v_{\theta})$

Part2.单位根

（以下，默认 $n$ 为 $2$ 的整数幂）

我们把单位圆上的某个复数

(c o s (2 k π / n), s i n (2 k π / n))

$(cos(2k\pi/n),sin(2k\pi/n))$ 称作单位根

ω_{n}^{k}

$\omega_n^k$ 。换句话说，就是极角表示法下的复数

(1, 2 k π / n)

$(1,2k\pi/n)$ 。当然，因为是在单位圆上，我们可以只关注它的辐角

(2 k π / n)

$(2k\pi/n)$ 。

如图。

![8次单位根](https://cdn.luogu.com.cn/upload/image_hosting/otcevc5c.png)

结论1：

ω_{n}^{k} = ω_{n}^{k + n}

$\omega_n^k=\omega_n^{k+n}$ 。

类比正负角。可以感性理解一下，就等于绕了一整圈，度数

> 2 π

$>2\pi$ 的角。有

ω_{n}^{k} = 2 k π / n = 2 k π / n + 2 π = 2 k π / n + 2 n π / n = 2 (k + n) π / n = ω_{n}^{k + n}

$\omega_n^k=2k\pi/n=2k\pi/n+2\pi=2k\pi/n+2n\pi/n=2(k+n)\pi/n=\omega_n^{k+n}$

结论2：

ω_{n}^{k} = - ω_{n}^{k + n / 2}

$\omega_n^k=-\omega_n^{k+n/2}$

等于绕了半圈的角。

结论3：当 $k$ 为偶数时， $\omega_n^k=\omega_{n/2}^{k/2}$ 。这是最重要的结论。

有

ω_{n}^{k} = 2 k π / n = 2 (k / 2) π / (n / 2) = ω_{n / 2}^{k / 2}

$\omega_n^k=2k\pi/n=2(k/2)\pi/(n/2)=\omega_{n/2}^{k/2}$

Part3.DFT

显然，

n + 1

$n+1$ 个点可以唯一确定一个

n

$n$ 次多项式。

而我们如果知道

f

$f$ 和

g

$g$ 上各

n + m + 1

$n+m+1$ 个点，就可以在

O (n)

$O(n)$ 时间內把它们乘在一起。对于每个点

(x, f (x))

$(x,f(x))$ 与

(x, g (x))

$(x,g(x))$ ，新点即为

(x, f (x) g (x))

$(x,f(x)g(x))$ 。

DFT解决的就是将多项式

f

$f$ 由系数表达转为点值表达。

开始推式子：

f (x)

$f(x)$

= a_{0} + a_{1} x + a_{2} x^{2} + . . . + a_{n} x^{n}

$=a_0+a_1x+a_2x^2+...+a_nx^n$

= (a_{0} + a_{2} x^{2} + . . . + a_{n} x^{n}) + (a_{1} x + a_{3} x^{3} + . . . + a_{n - 1} x^{n - 1})

$=(a_0+a_2x^2+...+a_nx^n)+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})$

= (a_{0} + a_{2} x^{2} + . . . + a_{n} x^{n}) + x (a_{1} + a_{3} x^{2} + . . . + a_{n - 1} x^{n - 2})

$=(a_0+a_2x^2+...+a_nx^n)+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})$

令

f 1 (x) = a_{0} + a_{2} x + . . . + a_{n} x^{n / 2}

$f1(x)=a_0+a_2x+...+a_nx^{n/2}$ ，

f 2 (x) = a_{1} + a_{3} x + . . . + a_{n - 1} x^{n / 2}

$f2(x)=a_1+a_3x+...+a_{n-1}x^{n/2}$

则

f (x) = f 1 (x^{2}) + x f 2 (x^{2})

$f(x)=f1(x^2)+xf2(x^2)$ 。

令

x = ω_{n}^{k} (k < n / 2)

$x=\omega_n^k(k<n/2)$ 。

则

f (ω_{n}^{k})

$f(\omega_n^k)$

= f 1 ((ω_{n}^{k})^{2}) + ω_{n}^{k} f 2 ((ω_{n}^{k})^{2})

$=f1((\omega_n^k)^2)+\omega_n^kf2((\omega_n^k)^2)$

= f 1 (ω_{n}^{2 k}) + ω_{n}^{k} f 2 (ω_{n}^{2 k})

$=f1(\omega_n^{2k})+\omega_n^kf2(\omega_n^{2k})$

= f 1 (ω_{n / 2}^{k}) + ω_{n}^{k} f 2 (ω_{n / 2}^{k})

$=f1(\omega_{n/2}^k)+\omega_n^kf2(\omega_{n/2}^k)$

再令

x = ω_{n}^{k + n / 2} (k < n / 2)

$x=\omega_n^{k+n/2}(k<n/2)$ 。

则

f (ω_{n}^{k + n / 2})

$f(\omega_n^{k+n/2})$

= f 1 ((ω_{n}^{k + n / 2})^{2}) + ω_{n}^{k + n / 2} f 2 ((ω_{n}^{k + n / 2})^{2})

$=f1((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}f2((\omega_n^{k+n/2})^2)$

= f 1 (ω_{n}^{2 k + n}) + ω_{n}^{k + n / 2} f 2 (ω_{n}^{2 k + n})

$=f1(\omega_n^{2k+n})+\omega_n^{k+n/2}f2(\omega_n^{2k+n})$

= f 1 (ω_{n}^{2 k}) + ω_{n}^{k + n / 2} f 2 (ω_{n}^{2 k})

$=f1(\omega_n^{2k})+\omega_n^{k+n/2}f2(\omega_n^{2k})$

= f 1 (ω_{n / 2}^{k}) - ω_{n}^{k} f 2 (ω_{n / 2}^{k})

$=f1(\omega_{n/2}^k)-\omega_n^kf2(\omega_{n/2}^k)$

注意到什么了吗？

f (ω_{n}^{k})

$f(\omega_n^k)$ 与

f (ω_{n}^{k + n / 2})

$f(\omega_n^{k+n/2})$ ，结果只有最后一项的符号不同！！！！！！**

然后就可以分治了。只要知道

f 1 (ω_{n / 2}^{k})

$f1(\omega_{n/2}^k)$ 与

f 2 (ω_{n / 2}^{k})

$f2(\omega_{n/2}^k)$ 的值，就可以直接得出

f (ω_{n}^{k})

$f(\omega_n^k)$ 与

f (ω_{n}^{k + n / 2})

$f(\omega_n^{k+n/2})$ ！！！

因此我们可以令

f 1 (ω_{n / 2}^{k})

$f1(\omega_{n/2}^k)$ 成为新的

f (x)

$f(x)$ ，进一步递推。

Part4.IDFT

只需要将DFT中的

ω_{n}^{k}

$\omega_n^k$ 改为

ω_{n}^{- k}

$\omega_n^{-k}$ 即可。别忘了最后记得除上数组长度！！！

Part5.代码:

#include<bits/stdc++.h>
using namespace std;
const int MAXN=4000000;
int n,m,lim=1,t,res[MAXN];
const double pi=acos(-1);
struct cp{
 double x,y;
 cp(double u=0.0,double v=0.0){
  x=u,y=v;
 }
 friend cp operator +(const cp &lv,const cp &rv){
  return cp(lv.x+rv.x,lv.y+rv.y);
 }
 friend cp operator -(const cp &lv,const cp &rv){
  return cp(lv.x-rv.x,lv.y-rv.y);
 }
 friend cp operator *(const cp &lv,const cp &rv){
  return cp(lv.x*rv.x-lv.y*rv.y,lv.x*rv.y+lv.y*rv.x);
 }
}f[MAXN],g[MAXN];
void FFT(cp *a,int sz,int inv){
 if(sz==1)return;
 int md=sz>>1;
 static cp b[MAXN];
 for(int i=0;i<md;i++)b[i]=a[i<<1],b[i+md]=a[(i<<1)+1];
 for(int i=0;i<sz;i++)a[i]=b[i];
 FFT(a,md,inv),FFT(a+md,md,inv);
 for(int i=0;i<md;i++){
  cp x=cp(cos(2*pi*i/sz),inv*sin(2*pi*i/sz));
  b[i]=a[i]+x*a[i+md],b[i+md]=a[i]-x*a[i+md];
 }
 for(int i=0;i<sz;i++)a[i]=b[i];
}
int read(){
 char c=getchar();
 int x=0;
 while(c>'9'||c<'0')c=getchar();
 while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
 return x;
}
int main(){
 n=read(),m=read();
 for(int i=0,t;i<=n;i++)f[i].x=read();
 for(int i=0,t;i<=m;i++)g[i].x=read();
 while(lim<=n+m)lim<<=1;
 FFT(f,lim,1),FFT(g,lim,1);
 for(int i=0;i<lim;i++)f[i]=f[i]*g[i];
 FFT(f,lim,-1);
 for(int i=0;i<lim;i++)res[i]=(int)(f[i].x/lim+0.5);
 for(int i=0;i<=n+m;i++)printf("%d ",res[i]);
 return 0;
}

Part6.非递归

在FFT时，一个位置上的数的最终位置，是把其位置的二进制表达翻转后的新位置。

即有:

(7)_{10} = (111)_{2} \to (111)_{2} = (7)_{10}

$(7)_{10}=(111)_2 \rightarrow(111)_2=(7)_{10}$

(13)_{10} = (1101)_{2} \to (1011)_{2} = (11)_{10}

$(13)_{10}=(1101)_2\rightarrow(1011)_2=(11)_{10}$

(22)_{10} = (10110)_{2} \to (01101)_{2} = (13)_{10}

$(22)_{10}=(10110)_2\rightarrow(01101)_2=(13)_{10}$

显然，当长度不同时，一个数可能翻转到不同的位置。

而它的终点位置，可以如此递推：

for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));

其中，

l i m

$lim$ 是长度，

l g = l o g_{2} (l i m)

$lg=log_2(lim)$ 。

Part7.最终非递归代码：

#include<bits/stdc++.h>
using namespace std;
const int MAXN=4000000;
int n,m,lim=1,t,res[MAXN],lg,rev[MAXN];
const double pi=acos(-1);
struct cp{
 double x,y;
 cp(double u=0.0,double v=0.0){
  x=u,y=v;
 }
 friend cp operator +(const cp &lv,const cp &rv){
  return cp(lv.x+rv.x,lv.y+rv.y);
 }
 friend cp operator -(const cp &lv,const cp &rv){
  return cp(lv.x-rv.x,lv.y-rv.y);
 }
 friend cp operator *(const cp &lv,const cp &rv){
  return cp(lv.x*rv.x-lv.y*rv.y,lv.x*rv.y+lv.y*rv.x);
 }
}f[MAXN],g[MAXN];
int read(){
 char c=getchar();
 int x=0;
 while(c>'9'||c<'0')c=getchar();
 while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
 return x;
}
void FFT(cp *a,int tp){
 for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
 for(int md=1;md<lim;md<<=1){
  cp rt=cp(cos(pi/md),tp*sin(pi/md));
  for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
   cp w=cp(1,0);
   for(int i=0;i<md;i++,w=w*rt){
    cp x=a[pos+i],y=w*a[pos+md+i];
    a[pos+i]=x+y;
    a[pos+md+i]=x-y;
   }
  }
 }
}
int main(){
 n=read(),m=read();
 while(lim<=(n+m))lim<<=1,lg++;
 for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
 for(int i=0;i<=n;i++)f[i].x=read();
 for(int i=0;i<=m;i++)g[i].x=read();
 FFT(f,1),FFT(g,1);
 for(int i=0;i<lim;i++)f[i]=f[i]*g[i];
 FFT(f,-1);
 for(int i=0;i<lim;i++)res[i]=(int)(f[i].x/lim+0.5);
 for(int i=0;i<=n+m;i++)printf("%d ",res[i]);
 return 0;
}

完结撒傅里叶~~~

![](https://cdn.luogu.com.cn/upload/image_hosting/4l5cvdy8.png)

posted @ 2020-04-22 21:48 Troverld 阅读(164) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型，支持深度思考和联网搜索！
· 基于 Docker 搭建 FRP 内网穿透开源项目（很简单哒）
· ollama系列1：轻松3步本地部署deepseek，普通电脑可用
· 按钮权限的设计及实现
· 【杂谈】分布式事务——高大上的无用知识？

Troverld

博客渲染出问题了那是博客的锅，本人什么都不会做的。

FFT感性瞎扯

FFT可以干啥？

Part1.复数

Part2.单位根

Part3.DFT

Part4.IDFT

Part5.代码:

Part6.非递归

Part7.最终非递归代码：

完结撒傅里叶~~~

公告

搜索

常用链接

随笔分类

随笔档案

阅读排行榜

评论排行榜

推荐排行榜