2. 程式人生 > >快速傅裏葉變換(FFT)

快速傅裏葉變換(FFT)

阿新 • • 發佈：2019-02-02
tdi 滿足 2.0 ++ n-2 amp 會有 ctime ast

扯

去北京學習的時候才系統的學習了一下卷積，當時整理了這個筆記的大部分。後來就一直放著忘了寫完。直到今天都臘月二十八了，才想起來還有個FFT的筆記沒整完呢。整理完這個我就假裝今年的任務全都over了吧。
更改了一些以前不大正確的地方，又添加了一些推導，證明實在不會。
有一些公式，但個人覺得還是比較好理解。可能還會有錯誤，希望大佬友情指出。
最後，祝各位看官新年快樂.
回家過寒假去咯(雖然就\(4\)天\(qwq\))

多項式

一個次數界為\(n\)的多項式\(A(x) = \sum_{i = 0}^{n - 1}a_ix^i\)
次數界為\(n\)說明這個多項式最高次項的指數小於\(n\)

系數表示法

我們可以只保留多項式每項的系數來描述這個多項式。即\(a_i\)表示指數為\(i\)的項的系數

運算

多項式之間的加運算，只要將對應項的系數相加即可。
即若\(C(x) = A(x) + B(x)\)那麽有\(C(x) = \sum_{i = 0}^{n - 1}{(a_i+b_i)x^i}\)
多項式之間的乘法。
若\(C(x) = A(x) + B(x)\)那麽\(C(x) = \sum_{i = 0}^{n - 1}{\sum_{j = 0}^ia_jb_{i - j}x^i}\)
其中\(C(x) = \sum_{i = 0}^{n - 1}{\sum_{j = 0}^ia_jb_{i - j}x^i}\)

叫做卷積
可以看出暴力計算卷積的復雜度為\(O(n^2)\)

點值表示法

我們可以選\(n\)個值代入多項式從而得出\(n\)個二元組\((x_i,y_i)\)。可以證明，通過這\(n\)個點對可以確定這個多項式。
還原這個多項式用高斯消元即可。
點值表示法的優點是進行多項式的加法和乘法只要將\(y_i\)相加或者相乘即可。

求值與插值

考慮如何可以比較快速的計算卷積。因為點值表示計算乘法比較方便。所以我們可以把系數表示法轉換為點值表示法。然後進行乘法。然後在轉換回來。就可以計算卷積了。
把系數表示法轉化為點值表示法稱為求值，把點值表示法轉換為系數表示法稱為插值。

拉格朗日插值公式

\(A(x) = \sum_{i = 0}^{n-1}y_i\frac{\prod_{j\ne i}{(x - x_j)}}{\prod_{j \ne i}{(x_i - x_j)}}\)


利用拉格朗日插值公式可以在給定\(n\)個點值表示法的情況下\(O(n^2)\)時間內計算出原二項式在某個位置的值

離散傅裏葉變換(DFT)

離散傅裏葉變換實現了兩種表示法之間的轉換，復雜度是\(O(n^2)\)

復數

復數是形如\(a+bi\)的數，其中\(a,b\)為實數。\(a\)稱為實部，\(b\)稱為虛部。\(i\)為虛數單位，定義\(i^2 = -1\)
復數是數的範圍在實數上的擴展
將復數\(a+bi\)表示在數軸上，是一個\((a,b)\)的向量.

復數的運算

復數相加將\(a\)和\(b\)分別相加即可。
復數相乘，其實將i看作未知量相乘就行了。即\((a+bi)*(x+yi) = (ax - by) + (ay + bx)i\)

單位圓

在數軸上以原點為圓心半徑為\(1\)的圓。

單位根

\(n\)次單位根滿足\(\omega^n = 1\)。\(\omega 就是n次單位根\)
表現在單位圓上就是單位圓的\(n\)等分點
我們用\(\omega_n\)表示主\(n\)次單位根也就是從\((1,0)\)開始數第一個n次單位根
那麽可以得到其他的\(n\)次單位根分別為\(\omega_n^0,\omega_n^2,\omega_n^3,.....\omega_n^{n - 1}\)
性質:
\(\omega_n^n = 1\)當n為偶數時,\(\omega_n^{n/2} = -1\)
\(\sum_{i=0}^{n-1}\omega_n^i=0\)

歐拉公式

\(e^{ix} = cos(x)+i * sin(x)\)
\(\omega_n=e^{2\pi i/n} = cos(\frac{2\pi}{n})+i*sin(\frac{2\pi}{n})\)
用來求主\(n\)次單位根

轉換

有了上面這些鋪墊，就可以進行轉換了。
將\(n\)個\(n\)次單位根帶入原來的多項式，就可以得到\(n\)個點值表示法。
那麽怎麽將點值表示法轉換為系數表示法呢？
這就是用\(n\)次單位根的原因
我們將得到\(n\)個\(y\)作為另一個多項式\(B(x)\)的系數
然後取單位根倒數帶入就能得到原來的多項式\(A(x)\)啦
這樣就成功的實現了點值表示法與系數表示法之間的轉換

快速傅裏葉變換(FFT)

用\(DFT\)計算卷積的復雜度是\(O(n^2)\)的。利用快速傅裏葉變換可以優化到\(O(nlogn)\)

分治

考慮一個次數界為\(n\)的多項式\(A(x)\)求他的\(DFT\)
這裏默認\(n\)為偶數，如果\(n\)不是偶數，那麽只要加上一個系數為\(0\)的項即可將他填充為偶數
然後我們根據下標的奇偶性將他分為兩部分。
\[A_{(even)}(x) = a_0 + a_2x + a_4x^2....+ a_{n-2}x^{\frac{n}{2} - 1}\]
\[A_{(odd)}(x) = a_1 + a_3x+a_4x^2....+a_{n-1}x^{\frac{n}{2}-1}\]
然後就能得到\(A(x)=A_{(even)}(x^2) + xA_{(odd)}(x^2)\)
然後發現\(A_{(even)}(x)\)和\(A_{(odd)}(x)\)是可以用同樣的方法遞歸計算的。所以就可以進行分治然後遞歸運算。
那麽我們需要帶入這\(n\)個單位根，應該怎麽計算呢。
下面進行推導
假如現在要把\(\omega_{n}^i\)代入
依據上面的式子
\[A(\omega_n^i)=A(\omega_n^{2i}) + \omega_n^{i}A(\omega_n^{2i})\]
\[=A(\omega_{\frac{n}{2}}^i)+\omega_n^iA(\omega_{\frac{n}{2}}^i)\]
然後考慮\(A(\omega_n^{i + \frac{n}{2}})\)
\[A(\omega_n^{i + \frac{n}{2}}) = A(\omega_{n}^{2i+n}) + \omega_n^{i + \frac{n}{2}} A(\omega_{n}^{2i+n})\]
\[=A(\omega_{n}^{2i}\omega_n^n) + \omega_n^{i} \omega_{n}^{ \frac{n}{2}} A(\omega_{n}^{2i}\omega_n^n)\]
\[=A(\omega_{\frac{n}{2}}^i) - \omega_n^{i}A(\omega_{\frac{n}{2}}^i)\]
然後就可以遞歸的愉快的計算出將這\(n\)個單位根代入所得到的值

優化

遞歸實現\(FFT\)非常慢。所以要用非遞歸的方法來實現
我們將遞歸實現的過程寫下來
0 1 2 3 4 5
0 2 4|1 3 5
0 4|2|1 5|3
0|4|2|1|5|3
然後看出每個數字遞歸到最後的為是為當前位置二進制表示法翻轉之後的數字。
然後就可以先把每個數字都放到最後的位置上面，然後再合並就行了。

代碼

遞歸

/*
* @Author: wxyww
* @Date:   2019-02-01 21:09:51
* @Last Modified time: 2019-02-01 21:47:34
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cmath>
#include<ctime>
#include<cmath>
#include<bitset>
using namespace std;
typedef long long ll;
const int N = 3000100;
const double pi = acos(-1.0);
ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
struct complex {
    double x,y;
    complex() {x = y = 0;}
    complex(double xx,double yy) {
        x = xx;y = yy;
    }
}A[N],B[N];
complex operator * (complex a,complex b) {
    return complex(a.x * b.x - a.y * b.y,a.x * b.y + a.y * b.x);//復數運算
}
complex operator + (complex a,complex b) {
    return complex(a.x + b.x,a.y + b.y);
}
complex operator - (complex a,complex b) {
    return complex(a.x - b.x,a.y - b.y);
}

void FFT(complex *a,int n,int ty) {
    if(n == 1) return;
    complex a1[n >> 1],a2[n >> 1];
    for(int i = 0;i <= n;i += 2) {
        a1[i >> 1] = a[i];a2[i >> 1] = a[i + 1];//按奇偶分開
    }
    FFT(a1,n >> 1,ty);FFT(a2,n >> 1,ty);//遞歸
    complex w1 = complex(cos(2.0 * pi / n),ty * sin(2.0 * pi / n));//主n次單位根
    complex w = complex(1.0,0.0);//當前的n次單位根
    int k = n >> 1;
    for(int i = 0;i < k;++i) {//根據式子計算
        complex t = w * a2[i];
        a[i + k] = a1[i] - t;
        a[i] = a1[i] + t;
        w = w *  w1;
    }

}
int main() {
    int n = read(),m = read();
    for(int i = 0;i <= n;++i) A[i].x = read();
    for(int i = 0;i <= m;++i) B[i].x = read();
    int tot = 1;
    while(tot <= n + m) tot <<= 1;
    FFT(A,tot,1);
    FFT(B,tot,1);
    for(int i = 0;i <= tot;++i) A[i] = A[i] * B[i];
    FFT(A,tot,-1);
    for(int i = 0;i <= n + m;++i) 
        printf("%d ",(int)(A[i].x / tot + 0.5));
    return 0;
}

非遞歸


/*
* @Author: wxyww
* @Date:   2019-02-01 21:50:51
* @Last Modified time: 2019-02-01 22:03:10
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cmath>
#include<ctime>
#include<cmath>
#include<bitset>
using namespace std;
typedef long long ll;
const int N = 3000100;
const double pi = acos(-1.0);
ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
struct complex {
    double x,y;
    complex() {x = y = 0;}
    complex(double xx,double yy) {
        x = xx;y = yy;
    }
}A[N],B[N];
complex operator * (complex a,complex b) {
    return complex(a.x * b.x - a.y * b.y,a.x * b.y + a.y * b.x);//復數運算
}
complex operator + (complex a,complex b) {
    return complex(a.x + b.x,a.y + b.y);
}
complex operator - (complex a,complex b) {
    return complex(a.x - b.x,a.y - b.y);
}

void FFT(complex *a,int n,int ty) {
    for(int i = 0,j = 0;i < n;++i) {//找到最終對應的位置
        if(i < j) swap(a[i],a[j]);
        for(int k = n >> 1;(j ^= k) < k;k >>= 1);
    }
    for(int m = 2;m <= n;m <<= 1) {
        complex w1 = complex(cos(2*pi/m),ty * sin(2 * pi / m));
        for(int i = 0;i < n;i += m) {
            complex w = complex(1,0);
            for(int k = 0;k < (m >> 1);++k) {
                complex t = w * a[i + k + (m >> 1)];
                complex u = a[i + k];
                a[i + k] = u + t;
                a[i + k + (m >> 1)] = u - t;
                w = w * w1;
            }
        }
    }
}
int main() {
    int n = read(),m = read();
    for(int i = 0;i <= n;++i) A[i].x = read();
    for(int i = 0;i <= m;++i) B[i].x = read();
    int tot = 1;
    while(tot <= n + m) tot <<= 1;
    FFT(A,tot,1);
    FFT(B,tot,1);
    for(int i = 0;i <= tot;++i) A[i] = A[i] * B[i];
    FFT(A,tot,-1);
    for(int i = 0;i <= n + m;++i) 
        printf("%d ",(int)(A[i].x / tot + 0.5));
    return 0;
}

時間差異

最後來一份遞歸代碼與非遞歸代碼運行時間的比較
技術分享圖片

快速傅裏葉變換(FFT)

相關推薦

快速傅裏葉變換(FFT)模板

while com lib rev using amp namespace for lap //有多少項N取多至少8~10倍 #include <algorithm> #include <cmath> #include <complex&

快速傅裏葉變換FFT

span n) 需要 單位根 空間 思想 mon clas splay FFT是基於分治思想對DFT和IDFT的優化。 DFT: \(\mathbb{C} ^{N}\rightarrow \mathbb{C} ^{N}: \left( x_{1},\ldots ,x_{

[模板]快速傅裏葉變換 FFT

nbsp namespace str using space clu ble friend getchar #include <iostream> #include <cstdio> #include <cstring> #in

[模板] 快速傅裏葉變換/fft

快速 cos tcp swap int -a 快速傅裏葉變換 ons urn Code const int nmax=(int)3e6+50; const db pi=acos(-1.0); struct tcpx{db a,b;}c1[nmax],c2[nmax]; t

快速傅裏葉變換(FFT)

tdi 滿足 2.0 ++ n-2 amp 會有 ctime ast 扯 去北京學習的時候才系統的學習了一下卷積，當時整理了這個筆記的大部分。後來就一直放著忘了寫完。直到今天都臘月二十八了，才想起來還有個FFT的筆記沒整完呢。整理完這個我就假裝今年的任務全都over了吧。

【Learning】多項式乘法與快速傅裏葉變換（FFT）

alt 技術 cos 相同 es2017 define ostream 意思 呵呵 簡介： 　　FFT主要運用於快速卷積，其中一個例子就是如何將兩個多項式相乘，或者高精度乘高精度的操作。 　　顯然暴搞是$O(n^2)$的復雜度，然而FFT可以將其將為$O(n lg

OI中的快速傅裏葉變換（FFT）

ali 系列 我認 交流 cpp 實現 AS 需要 部分 快速傅裏葉變換（FFT） 前言 關於這篇文章 ? ? 非常高興能有機會來探討快速傅裏葉變換，也就是大家熟知的 \(FFT\) 在 \(OI\) 中的運用。以前了解過一次 \(FFT\) ，現在過了幾個月，數學和 \

洛谷P1919 A*B problem 快速傅裏葉變換模板 [FFT]

sca 快速 https -o ora stream ros for AR 　　題目傳送門 A*B problem 題目描述 給出兩個n位10進制整數x和y，你需要計算x*y。 輸入輸出格式 輸入格式： 第一行一個正整數n。 第二行描述一個位數為n的正整數

【FFT】快速傅裏葉變換

運算 lse ssi ive b+ mathjax ret 方向 strong 【FFT】快速傅裏葉變換 一、復數 1、定義 復數：設 $a$，$b$ 為實數，$i^{2}=−1$ ，形如 $a+bi$ 的數叫復數，其中 $i$ 被稱為虛數單位，復數域是目

【模板】快速傅裏葉變換（FFT）（BZOJ2179）

turn const 一個 algorithm har string using tchar out Description 給出兩個n位10進制整數x和y，你需要計算x*y。 Input 第一行一個正整數n。 第二行描述一個位數為n的正整數x。 第三行描述一個位數為

【bzoj2179】FFT快速傅裏葉變換（優化高精度乘法）

efi swa 快速傅裏葉變換 得到 urn lse can sin %s #include<bits/stdc++.h> using namespace std; #define pi acos(-1) typedef complex<double&g

FFT 快速傅裏葉變換 學習筆記

.net 消元 esp line tac ace 表示 etc cnblogs FFT 快速傅裏葉變換 前言 lmc,ikka,attack等眾多大佬都沒教會的我終於要自己填坑了。 又是機房裏最後一個學fft的人 早背過圓周率50位填坑了 用處 多項式乘法 卷積 \(g(x

快速傅裏葉變換

idt ima .com com ges 分享 logs jpg cnblogs 快速傅裏葉變換

Codeforces 827E Rusty String - 快速傅裏葉變換 - 暴力

後綴 graph pro 問題 luci sid alt image .com Grigory loves strings. Recently he found a metal strip on a loft. The strip had lengt

51Nod 快速傅裏葉變換題集選刷

哈希 getc href 圖片 ase 推出 hid 階乘 序列 打開51Nod全部問題頁面，在右邊題目分類中找到快速傅裏葉變換，然後按分值排序，就是本文的題目順序。 1.大數乘法問題 這個……板子就算了吧。 2.美妙的序列問題

真實感海洋的繪制（二）：使用快速傅裏葉變換加速波形計算

image 完全 eps gpu spa 可能 src 重寫 時間 真實感海洋的繪制（二）：使用快速傅裏葉變換加速波形計算 其實上一篇博文所寫的\(H(\vec{x},t)?\)，就是二維傅裏葉變換的求和式，之前的暴力計算法屬於二維的離散傅裏葉變換（Discrete Fou

【模板】快速傅裏葉變換

wap body problem rev pan pos bit urn 傅裏葉變換 uoj34 1 #include<bits/stdc++.h> 2 #define db double 3 using namespace std; 4 const

caioj1450:【快速傅裏葉變換】大整數乘法

rose clas scan name 代碼 printf 答案 r+ 傅裏葉變換 【傳送門：caioj1450】 簡要題意： 　　給出兩個超級大的整數，求出a*b 題解： 　　Rose_max出的一道FFT例題，卡掉高精度 = = 　　只要把a和b的每一

快速傅立葉變換FFT模板

return namespace double names http ++ main swap pre 遞歸版 UOJ34多項式乘法 //容易暴棧，但是很好理解 #include <cmath> #include <iostream> #includ

快速傅立葉變換FFT的學習筆記一：C語言程式碼的簡單實現

快速傅立葉變換FFT的學習筆記一：C語言程式碼的簡單實現 fft.c #include "math.h" #include "fft.h" void conjugate_complex(int n,complex in[],complex out[]) { int i = 0