整理的演算法模板合集： ACM模板
目錄
- 多項式
- 復數
- 復數的單位根 / 單位向量
- 離散傅里葉變換（DTF）
- 離散傅里葉逆變換（IDTF）
- FFT演算法整體流程
- 代碼實作
- 遞回版
- 迭代版
- “三步變兩步”優化
- P1919 【模板】FFT快速傅里葉變換
學習筆記（用ipad記錄的筆記）：
超級簡單的快速傅里葉變換！只要基礎夠扎實，順著推一遍沒有什么難以理解的，我學的整個程序沒有一點卡殼，真的很爽，整整寫了三個多小時，寫了滿滿6頁的筆記，主要是內容太多了，
首先是一些基礎概念：
多項式
復數
復數的單位根 / 單位向量
離散傅里葉變換（DTF）
離散傅里葉逆變換（IDTF）
FFT演算法整體流程
代碼實作
遞回版
我們把上述的思路實作用遞回可以很形象地實作FFT函式，
這里要注意一下，盡管C++自帶有復數庫，但是建議手寫一下，代碼不長，手寫常數小，
struct Complex
{
double x, y;
Complex (double x = 0, double y = 0) : x(x), y(y) { }
}A[N], B[N];
Complex operator * (Complex J, Complex Q) {
//模長相乘，幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator - (Complex J, Complex Q) {
return Complex(J.x - Q.x, J.y - Q.y);
}
Complex operator + (Complex J, Complex Q) {
return Complex(J.x + Q.x, J.y + Q.y);
}
void FFT(int limit, Complex *a, int type) {
if (limit == 1) return ; //只有一個常數項
Complex a1[limit >> 1], a2[limit >> 1];
for (int i = 0; i <= limit; i += 2) //根據下標的奇偶性分類
a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
FFT(limit >> 1, a1, type);
FFT(limit >> 1, a2, type);
Complex Wn = Complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = Complex(1, 0);
//Wn為單位根，w表示冪
for (int i = 0; i < (limit >> 1); i++, w = w * Wn) //這里的w相當于公式中的k
a[i] = a1[i] + w * a2[i],
a[i + (limit >> 1)] = a1[i] - w * a2[i]; //利用單位根的性質，O(1)得到另一部分
}
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 limit = 1; while (limit <= N + M) limit <<= 1;
FFT(limit, a, 1);
FFT(limit, b, 1);
//后面的1表示要進行的變換是什么型別
//1表示從系數變為點值
//-1表示從點值變為系數
//至于為什么這樣是對的，可以參考一下c向量的推導程序，
for (int i = 0; i <= limit; i++)
a[i] = a[i] * b[i];
FFT(limit, a, -1);
for (int i = 0; i <= N + M; i++) printf("%d ", (int)(a[i].x / limit + 0.5)); //按照我們推倒的公式，這里還要除以n
return 0;
}
但是這樣寫常數太大了，還需要一個很大的陣列存，實測會T飛
迭代版
所以我們引入一種迭代版，利用蝴蝶變換巧妙地實作FFT，使其真正成為 O ( n l o g n ) O(nlogn) O(nlogn)的高效演算法，（下面的幾張截圖摘自這里）
int bit = log2(N) - 1;
for (int i = 0; i < N; ++i)
{
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit);
if (i < rev[i]) // 這里要防止防止重復交換
swap(a[i], a[rev[i]]);
}
然后下面就是FFT的迭代模板：
P3803 【模板】多項式乘法（FFT）
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int MAXN = 1e7 + 10;
inline int read() {
char c = getchar();
int x = 0, f = 1;
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;
}
const double Pi = acos(-1.0);
struct Complex {
double x, y;
Complex(double xx = 0, double yy = 0) {
x = xx, y = yy;
}
} a[MAXN], b[MAXN];
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);
}
Complex operator * (Complex a, Complex b) {
return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); //不懂的看復數的運算那部分
}
int N, M;
int l, r[MAXN];
int limit = 1;
void FFT(Complex *A, int type) {
for (int i = 0; i < limit; i++)
if (i < r[i])
swap(A[i], A[r[i]]); //求出要迭代的序列
for (int mid = 1; mid < limit; mid <<= 1) { //待合并區間的長度的一半
Complex Wn(cos(Pi / mid), type * sin(Pi / mid)); //單位根
for (int R = mid << 1, j = 0; j < limit; j += R) { //R是區間的長度，j表示前已經到哪個位置了
Complex w(1, 0); //冪
for (int k = 0; k < mid; k++, w = w * Wn) { //列舉左半部分
Complex x = A[j + k], y = w * A[j + mid + k]; //蝴蝶效應
A[j + k] = x + y;
A[j + mid + k] = x - y;
}
}
}
}
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();
while (limit <= N + M)
limit <<= 1, l++;
for (int i = 0; i < limit; i++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)) ;
// 在原序列中 i 與 i/2 的關系是 ： i可以看做是i/2的二進制上的每一位左移一位得來
// 那么在反轉后的陣列中就需要右移一位，同時特殊處理一下奇數
FFT(a, 1);
FFT(b, 1);
for (int i = 0; i <= limit; i++)
a[i] = a[i] * b[i];
FFT(a, -1);
for (int i = 0; i <= N + M; i++)
printf("%d ", (int)(a[i].x / limit + 0.5));
return 0;
}
“三步變兩步”優化
設 P P P和 Q Q Q 是 實多項式， F = P + Q i F = P+Qi F=P+Qi ，則 F 2 = P 2 ? Q 2 + 2 P Q i F^2 = P^2-Q^2+2PQi F2=P2?Q2+2PQi ，注意到我們要求的 P Q PQ PQ正是 F F F 虛部的一半，這樣只需要兩次FFT就可以求出結果，
for (int i = 0; i <= n; ++i)
A[i] = read();
for (int i = 0; i <= m; ++i)
B[i] = read();
for (int i = 0; i <= max(n, m); ++i)
F[i] = comp(A[i], B[i]);
fft(F, N);
for (int i = 0; i < N; ++i)
F[i] = F[i] * F[i];
fft(F, N, -1);
for (int i = 0; i <= n + m; ++i)
printf("%d ", int(F[i].imag() / (N * 2) + 0.1));
P1919 【模板】FFT快速傅里葉變換
題目大意：給你兩個整數，a和b，求 a × b a\times b a×b， 其中 1 ≤ a , b ≤ 1 0 1000000 1≤a,b≤10^{1000000} 1≤a,b≤101000000，
對于每一個 n 位的十進制數，我們都可以看做一個 n-1 次多項式 A，滿足
A ( x ) = a 0 + a 1 × 10 + a 2 × 1 0 2 ? + a n ? 1 × 1 0 n ? 1 A(x) =a_0+a_1 \times 10+a_2\times10^2 \cdots +a_{n-1}\times10^{n-1} A(x)=a0?+a1?×10+a2?×102?+an?1?×10n?1
然后對于這兩個非常大的數，我們就可以用FFT快速求解答案了，
-
不要忘了進位！
-
不要忘了要保證數位上的單調性！因為我們普通的FFT卷積時，高次項一定由低次項得到，放在這里也一樣,所以我們要倒序存盤，
#include <cstdio>
#include <cmath>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N = 5000007;
const double PI = acos(-1);
int n, m;
int res, ans[N];
int AA, BB;
int Lim = 1;//
int L;//二進制的位數
int R[N];
inline int read()
{
register int x = 0, f = 1;
register char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-')f = -1;ch =getchar();}
while(ch >= '0' && ch <= '9') {x = x * 1- + ch - '0';ch = getchar();}
return x * f;
}
struct Complex
{
double x, y;
Complex (double x = 0, double y = 0) : x(x), y(y) { }
}A[N], B[N];
Complex operator * (Complex J, Complex Q) {
//模長相乘，幅度相加
return Complex(J.x * Q.x - J.y * Q.y, J.x * Q.y + J.y * Q.x);
}
Complex operator - (Complex J, Complex Q) {
return Complex(J.x - Q.x, J.y - Q.y);
}
Complex operator + (Complex J, Complex Q) {
return Complex(J.x + Q.x, J.y + Q.y);
}
string s1, s2;
inline void FFT(Complex *J, double type)
{
for(int i = 0; i < Lim; ++ i) {
if(i < R[i]) swap(J[i], J[R[i]]);
//i小于R[i]時才交換，防止同一個元素交換兩次，回到它原來的位置，
}
//從底層往上合并
for(int mid = 1; mid < Lim; mid <<= 1) {//待合并區間長度的一半，最開始是兩個長度為1的序列合并,mid = 1;
Complex wn(cos(PI / mid), type * sin(PI / mid));//單位根w_n^i;
for(int len = mid << 1, pos = 0; pos < Lim; pos += len) {
//for(int pos = 0; pos < Lim; pos += (mid << 1)) {
//len是區間的長度，pos是當前的位置,也就是合并到了哪一位
Complex w(1, 0);//冪,一直乘，得到平方，三次方...
for(int k = 0; k < mid; ++ k, w = w * wn) {
//只掃左半部分，得到右半部分的答案,w 為 w_n^k
Complex x = J[pos + k];//左半部分
Complex y = w * J[pos + mid + k];//右半部分
J[pos + k] = x + y;//蝴蝶效應
J[pos + mid + k] = x - y;
}
}
}
}
int cnt_a, cnt_b;
int main()
{
cin >> s1 >> s2;
n = s1.length();
m = s2.length();
//相當于x=10 的多項式
//讀入的數的每一位看成多項式的一項，保存在復數的實部
for (int i = n - 1; i >= 0; -- i) A[cnt_a ++ ].x = s1[i] - 48;
for (int i = m - 1; i >= 0; -- i) B[cnt_b ++ ].x = s2[i] - 48;
while(Lim < n + m) Lim <<= 1, L ++ ;
for (int i = 0; i <= Lim; ++ i) {
//換成二進制序列
R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
// 在原序列中 i 與 i/2 的關系是 ： i可以看做是i/2的二進制上的每一位左移一位得來
// 那么在反轉后的陣列中就需要右移一位，同時特殊處理一下奇數
}
//FFT 把a的系數表示轉化為點值表示
FFT(A, 1);
//FFT 把b的系數表示轉化為點值表示
FFT(B, 1);
//兩個多項式的點值表示相乘
for (int i = 0; i <= Lim; ++ i) A[i] = A[i] * B[i];
//IFFT 把這個點值表示轉化為系數表示
FFT(A, -1);
int tot = 0;
//保存答案的每一位（注意進位）
for (int i = 0; i <= Lim; ++ i) {
//取實數四舍五入，此時虛數部分應當為0或由于浮點誤差接近0
ans[i] += (int) (A[i].x / Lim + 0.5);
if(ans[i] >= 10) {
ans[i + 1] += ans[i] / 10;
ans[i] %= 10;
Lim += (i == Lim);//進一位
}
}
//刪掉前導零
while(!ans[Lim] && Lim >= 1) Lim -- ;
Lim ++ ;
while(-- Lim >= 0) cout << ans[Lim];
puts("");
return 0;
}
參考資料
-
快速傅里葉變換(FFT)詳解
-
演算法學習筆記(32): 快速傅里葉變換
-
一小時學會快速傅里葉變換（Fast Fourier Transform）
-
FFT·快速傅立葉變換
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/240567.html
標籤：其他