FFT

前置知识

两种多项式表示法:

系数表示法: A(x)=n1i=0aixi A ( x ) = ∑ i = 0 n − 1 a i x i ,其中 {an1} { a n − 1 } 就是系数。

点值表示法:代入 n n 个不同的 x 得到对应的 A(x) A ( x ) 记为 (x,y) ( x , y ) ,这 n n 个点 (x0,y0),(x1,y1)(xn1,yn1) 唯一确定了 A(x) A ( x )

由系数表示法变为点值表示法称为求值,点值表示法变为系数表示法称为插值。

点值表示法好处:快速求 A(x)×B(x) A ( x ) × B ( x ) 的点值表达式——只需要将相同 x x y 乘起来。

单位复数根:

n n 个模为 1 的复数 n n 等分复平面,令幅角最小的是 ωn ( n n 次单位复数根),那么其他 n1 个就是 ω2n,ω3nωnn ω n 2 , ω n 3 ⋯ ω n n (两个复数相乘,模长相乘,幅角相加),画图yy一下可以发现 ω0n=ωnn=1 ω n 0 = ω n n = 1

因为是 n n 等分,所以 ωnk=(cos2kπn,sin2kπn) ,由此可见 ω2k2n=ωkn,ωk+n2n=ωkn ω 2 n 2 k = ω n k , ω n k + n 2 = − ω n k

快速傅里叶变换

因为点值表示法可以快速进行多项式乘法,所以我们想将系数表示法转换成点值表示法,再进行多项式相乘,最后将点值表示法转换成系数表示法。两个过程分别为离散傅里叶变换和傅里叶逆变换。

ps:接下来都默认 n n 2k (如果不够那么高位补 0 0 即可)。

离散傅里叶变换:

如果将 A(x) 的系数奇偶拆分,我们可以得到:

A0(x)=a0+a2x+a4x2++an2xn21A1(x)=a1+a3x+a5x2++an1xn21A(x)=A0(x2)+xA1(x2) A 0 ( x ) = a 0 + a 2 x + a 4 x 2 + ⋯ + a n − 2 x n 2 − 1 A 1 ( x ) = a 1 + a 3 x + a 5 x 2 + ⋯ + a n − 1 x n 2 − 1 A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 )

因为要代入 n n 个不同的 x ,我们会联想到 n n 次单位复数根,又因为如果 n 是偶数,则 (ωkn)2=ωkn2 ( ω n k ) 2 = ω n 2 k ,也就是说 (ωkn)2 ( ω n k ) 2 的取值只有 n2 n 2 种。于是我们惊奇的发现只要代入 n n n 次单位复数根, A0(x2) A 0 ( x 2 ) A1(x2) A 1 ( x 2 ) 就变成了与 A(x) A ( x ) 同样的问题(这就是为什么 n n 要是 2k )!

所以可以用 O(nlog2n) O ( n l o g 2 n ) 的复杂度来愉快的离散傅里叶变换。

傅里叶逆变换:

只需要用 ωkn ω n − k 再做一次离散傅里叶变换,证明可以参考Menci大佬的博客QAQ,我太弱了就只记了结论:

ak=1ni=0n1ai(ωkn)i a k = 1 n ∑ i = 0 n − 1 a i ( ω n − k ) i

实现

离散傅里叶变换的递归写法很显然,但是常数巨大啊。实际上有迭代写法,我们观察每个 i i <script type="math/tex" id="MathJax-Element-48">i</script> 递归到最后的位置:

0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7
二进制:
 0   1   2   3   4   5   6   7
000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111
 0   4   2   6   1   5   3   7

惊奇的发现竟然是二进制反转QAQ,这样我们就可以从底向上迭代处理了。

模板

UOJ34,烂大街的FFT模板题。

#include<cstdio>
#include<cctype>
#include<cmath>
#include<algorithm>
#define fr first
#define sc second
#define mp make_pair
using namespace std;
typedef long double DB;typedef pair<DB,DB> C;
const int maxn=262144;const DB pi=acos(-1);

int n,m,R[maxn+5];C a[maxn+5],b[maxn+5];

#define Eoln(x) ((x)==10||(x)==13||(x)==EOF)
inline char readc(){
    static char buf[100000],*l=buf,*r=buf;
    if (l==r) r=(l=buf)+fread(buf,1,100000,stdin);
    if (l==r) return EOF;return *l++;
}
inline int readi(int &x){
    int tot=0,f=1;char ch=readc(),lst='+';
    while (!isdigit(ch)) {if (ch==EOF) return EOF;lst=ch;ch=readc();}
    if (lst=='-') f=-f;
    while (isdigit(ch)) tot=(tot<<3)+(tot<<1)+ch-48,ch=readc();
    return x=tot*f,Eoln(ch);
}
inline int Rev(int x,int len){
    static int buf[31];for (int i=0;i<len;i++) buf[i]=x&1,x>>=1;
    for (int i=0;i<len;i++) x=x<<1|buf[i];return x;
}
C operator + (const C &a,const C &b) {return mp(a.fr+b.fr,a.sc+b.sc);}
C operator - (const C &a,const C &b) {return mp(a.fr-b.fr,a.sc-b.sc);}
C operator * (const C &a,const C &b) {return mp(a.fr*b.fr-a.sc*b.sc,a.fr*b.sc+a.sc*b.fr);}
inline void FFT(C *a,int n,int f){
    for (int i=0;i<n;i++) if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<n;k<<=1){
        C w=mp(1,0),wn=mp(cos(pi/k),sin(f*pi/k)),x,y;
        for (int i=0;i<n;i+=(k<<1),w=mp(1,0))
            for (int j=0;j<k;j++,w=w*wn)
                x=a[i+j],y=w*a[i+j+k],a[i+j]=x+y,a[i+j+k]=x-y;
    }
}
int main(){
    freopen("FFT.in","r",stdin);
    freopen("FFT.out","w",stdout);
    readi(n);readi(m);
    for (int i=0,x;i<=n;i++) readi(x),a[i]=mp(x,0);
    for (int i=0,x;i<=m;i++) readi(x),b[i]=mp(x,0);
    m+=n;int len=0;for (n=1;n<=m;n<<=1) len++;for (int i=0;i<n;i++) R[i]=Rev(i,len);
    FFT(a,n,1);FFT(b,n,1);for (int i=0;i<n;i++) a[i]=a[i]*b[i];
    FFT(a,n,-1);for (int i=0;i<=m;i++) printf("%d ",(int)(a[i].fr/n+0.1));
    return 0;
}
  • 0
    点赞
  • 0
    收藏
  • 打赏
    打赏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:编程工作室 设计师:CSDN官方博客 返回首页
评论 1

打赏作者

ZigZagK

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值