作者：中二攻子

鏈接：https://ac.nowcoder.com/discuss/175409

來源：牛客網

本文含NTT、MTT、拆系數FFT、共軛優化FFT、多項式求逆與ln

約定：

1、

表示一個普通的項數為 的冪次多項式， 是他的點值表示。

2、

代表單位根， 表示 次單位根。

3、

代表一個數列。

4、

表示原根。

多項式系列之零 底層知識：

多項式的表示：

多項式可以通過系數數列

表示， 是 的系數。

多項式可以通過點值表示，對于一個

次多項式，取 種不同的

取值帶入

，得到 個值，在取相同這 個數的意義下，可以唯一的表示這個多項式。

多項式乘法：

定義

，在系數表示之下相乘復雜度 ，在點值表示之下 ，復雜度 。

復數：

復數一般情況下可以表示成

的形式， 是實數， 。

復數的幅角：平面直角坐標系上點

所在的任意角。

復數的模長：

兩個復數相乘：

，復數相乘之后，模長等于原來兩個復數的模長的乘積，幅角的角度等于原來兩個幅角的和。

復數可以加減乘除，可以和實數一樣的帶入

。

：

在單位圓上從

開始平均取 個點，從 開始編號，分別是 。

畫圖觀察可得：

所代表的復數

所代表的復數

DFT&IDFT：

在科學的數學函數意義上DFT是講一個函數轉化成三角函數的加減乘除的形式，三角函數的系數是原函數系數與點值之間的變換規律。IDFT是DFT的逆變換。

：

1、什么是

：在 意義下 互不相同，即 可以張成整個 下的域。

2、

存在的條件： ， 是奇素數。

3、如何求

：把 進行質因數分解 ，如果對于任意的 ，總有 ，暴力枚舉即可。

CRT合并：

求解

由

，得

帶入二式，得

若

，用逆元直接除便可；否則通過 可求得 ，若無解則方程組無解。

最后

。

多項式全集之一 FFT：

什么是FFT：

FFT是利用DFT的特殊性質，把

帶入 從而 求一個系數多項式的點值表示，所以叫FDFT。

的具體應用：

1、可以方便的IDFT：

設

的系數是 ，在 的DFT下點值是 ， 的系數是 ，在 的DFT下點值是 。

當

時 ，否則根據等比數列求和公式得

由此可得：

，

綜上所述，對于點值取的

相反數做DFT再除以 可得到系數。

2、可以快速的DFT：

直接將

帶入多項式做DFT需要復雜度 ，我們利用 的性質優化：

把

按照奇偶分裂，

我們令

我們可以發現

現在我們把

帶入，令 ：

我們知道取

時， 的取值，就可以算出 ，而 的長度都為 的一半，所以可以遞歸計算。

非遞歸優化FFT：

1、優化原理：

畫圖可知，遞歸版FFT最底層結束狀態第

個位置的項是 二進制翻轉后的結果。我們可以 的得到最底層的結果，然后向上模擬回溯合并即可。

2、蝴蝶變換：

由上述式子：

可得在迭代時

與 都只與 有關，所以我們可以用臨時變量記錄下一層的兩個信息向上迭代。

共軛復數優化FFT：

1、優化原理：

（在DFT時）

令

那么

2、證明：

：

：為方便起見，我們用 代替

而在IDFT時，我們需要

數論優化FFT(NTT)：

與 的共性：

1、

和 都互不相同 。

2、

， 。

3、

，由于原根的互不相同， ， 。

4、

，

5、

，

因為有這些共性，所以

可以代替 。

喜聞樂見的模板：

FFT模板（共軛優化）

namespace FFT{const double pi = acos(-1);struct cp{double x, y;cp() {x = y = 0;}cp(double X,double Y) {x = X; y = Y; }cp conj() {return (cp) {x, -y};}}a[3000005], b[3000005], c[3000005], I(0, 1), d[3000005];cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}cp operator* (const cp &a, double b) {return (cp){a.x * b, a.y * b};}cp operator/ (const cp &a, double b) {return (cp){a.x / b, a.y / b};}struct p_l_e{int wz[3000005];void FFT(cp *a, int N, int op) {for(int i = 0; i < N; i++)if (i<wz[i]) swap(a[i],a[wz[i]]);for(int le = 2; le <= N; le <<= 1) {int mid = le >> 1;for(int i = 0; i < N; i += le) {cp x, y, w = (cp) {1, 0};cp wn = (cp){cos(op * pi / mid), sin(op * pi / mid)};for(int j = 0 ; j < mid; j++) {x = a[i + j]; y = a[i + j + mid] * w;a[i + j] = x + y;a[i + j + mid] = x - y;w = w * wn;}}}}void D_FFT(cp *a, cp *b, int N, int op){for(int i = 0; i < N; i++) d[i] = a[i] + I * b[i];FFT(d, N, op);d[N] = d[0];for(int i = 0; i < N; i++){a[i] = (d[i] + d[N - i].conj()) / 2;b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;}d[N] = cp(0, 0);}void mult(cp *a, cp *b, cp *c, int M){int N = 1, len = 0;while(N < M) N <<= 1, len++;for(int i = 0; i < N; i++)wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));D_FFT(a, b, N, 1);for(int i = 0; i < N; i++) c[i] = a[i] * b[i];FFT(c, N, -1);for(int i = 0; i < N; i++) c[i].x = c[i].x / N;}}PLE;int n, m;void main() {scanf("%d%d", &n, &m); n++; m++;for(int i = 0; i < n; i++) scanf("%lf", &a[i].x);for(int i = 0; i < m; i++) scanf("%lf", &b[i].x);PLE.mult(a, b, c, n + m - 1); for(int i = 0; i < n + m - 1; i++) printf("%d ", (int)round(c[i].x));return;} }

NTT模板：

namespace NTT{typedef long long LL;const int mod = 998244353;const int g = 3;LL a[3000005], b[3000005], c[3000005];int n, m;LL qpow(LL a, LL b){LL ans = 1;while(b){if(b & 1) ans = ans * a % mod;a = a * a % mod;b >>= 1;}return ans;}struct p_l_e{int wz[3000005];void NTT(LL *a, int N, int op) {for(int i = 0; i < N; i++) if(i < wz[i]) swap(a[i], a[wz[i]]);for(int le = 2; le <= N; le <<= 1) {int mid = le >> 1;LL wn = qpow(g, (mod - 1) / le);if(op == -1) wn = qpow(wn, mod - 2); for(int i = 0; i < N; i += le) {int w = 1, x, y;for(int j = 0; j < mid; j++) {x = a[i + j]; y = a[i + j + mid] * w % mod; a[i + j] = (x + y) % mod;a[i + j + mid] = (x - y + mod) % mod;w = w * wn % mod;}}}}void mult(LL *a, LL *b, LL *c, int M) {int N = 1, len = 0;while(N < M) N <<= 1, len++;for(int i = 0; i < N; i++) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));NTT(a, N, 1); NTT(b, N, 1);for(int i = 0; i < N; i++) c[i] = a[i] * b[i] % mod;NTT(c, N, -1);LL t = qpow(N, mod - 2);for(int i = 0; i < N; i++) c[i] = c[i] * t % mod;}}PLE;void main() {scanf("%d%d", &n, &m); n++; m++;for(int i = 0; i < n; i++) scanf("%lld", &a[i]);for(int i = 0; i < m; i++) scanf("%lld", &b[i]);PLE.mult(a, b, c, n + m - 1);for(int i = 0; i < n + m - 1; i++) printf("%lld ", c[i]);} }

多項式全集之二 任長任模的FFT:

三模NTT實現任模FFT：

1、為什么要用MTT：當

不是NTT模數或者多項式長度大于模數限制時，就要使用MTT。

2、MTT的使用原理：我們對初始多項式取模，那么如果在不取模卷積情況下，答案

不會超過 。我們取三個NTT模數 ，分別做多項式乘法，得到 分別 的答案，通過CRT合并可以得到 的答案，如果 那么就可以得到準確的答案，再對 取模即可。

3、CRT合并的小優化：

初始式子

把一式二式合并(LL范圍內)。

再次合并(不需要 快速乘)。

4、常用NTT模數：

以下模數的共同

拆系數FFT(CFFT)實現任模FFT：

1、實現原理：運用實數FFT不取模做乘法，然后取模回歸到整數。但是由于誤差較大（值域是

），我們令

把系數 ，對 交叉做四遍卷積，求出答案按系數貢獻取模加入。

2、可按合并DFT的方法優化DFT次數。

算法實現任長FFT：

當

不是 的冪次的時候，我們從式子入手：

令

喜聞樂見的模板：

三模NTT模板（注意：不可以MTT回來，因為系數會取模）

namespace MTT{typedef long long LL;int n, m;LL p, mod;const LL p1 = 998244353;const LL p2 = 1004535809;const LL p3 = 104857601;const int g = 3189;LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005];LL c3[300005], c1[300005], c2[300005];LL qpow(LL a, LL b, LL mod) {LL ans = 1;while(b) {if(b & 1) ans = ans * a % mod;a = a * a % mod;b >>= 1;}return ans;}const LL inv12 = qpow(p1, p2 - 2, p2);const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3);struct p_l_e{int wz[300005];void MTT(LL *a, int N, int op) {for(int i = 0; i < N; i++)if(i < wz[i]) swap(a[i], a[wz[i]]);for(int le = 2; le <= N; le <<= 1) {int mid = le >> 1;LL wn = qpow(g, (mod - 1) / le, mod);if(op == -1) wn = qpow(wn, mod - 2, mod);for(int i = 0; i < N ;i += le) {LL w = 1, x, y;for(int j = 0; j < mid; j++) {x = a[i + j];y = a[i + j + mid] * w % mod;a[i + j] = (x + y) % mod;a[i + j + mid] = (x - y + mod) % mod;w = w * wn % mod;}}}}void mult(LL *a, LL *b, LL *c, int M) {int N = 1, len = 0;while(N < M) N <<= 1, len++;for(int i = 0; i < N; i++)wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));MTT(a, N, 1); MTT(b, N, 1);for(int i = 0; i < N; i++) c[i] = a[i] * b[i] % mod;MTT(c, N, -1);LL t = qpow(N, mod - 2, mod);for(int i = 0; i < N; i++) c[i] = c[i] * t % mod;}}PLE;LL CRT(LL c1, LL c2, LL c3) {LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2));LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p;return y;}void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) {for(int i = 0; i < N; i++)c[i] = CRT(c1[i], c2[i], c3[i]);return;}void main() {scanf("%d%d%lld", &n, &m, &p); n++; m++;for(int i = 0; i < n; i++) scanf("%lld", &a[i]);for(int i = 0; i < m; i++) scanf("%lld", &b[i]);mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1);mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1);mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1);merge(c1, c2, c3, c, n + m - 1);for(int i = 0; i < n + m - 1; i++) printf("%lld ", (c[i] % p + p) % p);return;} }

拆系數FFT模板（注意：相同系數的兩項可以合并一起IDFT。采用共軛優化法，只進行四次DFT）

namespace CFFT{typedef long long LL;int n, m, p ,sqrp; int a[300005], b[300005];const long double pi = acos(-1);struct cp{long double x, y;cp() {x = y = 0;}cp(long double X,long double Y) {x = X; y = Y; }cp conj() {return (cp) {x, -y};}}ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005];cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};}cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};} struct p_l_e{int wz[300005];void FFT(cp *a, int N, int op){for(int i = 0; i < N; i++)if(i < wz[i]) swap(a[i], a[wz[i]]);for(int le = 2; le <= N; le <<= 1){int mid = le >> 1;cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)};for(int i = 0; i < N; i += le){w = (cp){1, 0};for(int j = 0; j < mid; j++){x = a[i + j];y = a[i + j + mid] * w;a[i + j] = x + y;a[i + j + mid] = x - y;w = w * wn;}}} }void D_FFT(cp *a, cp *b, int N, int op){for(int i = 0; i < N; i++) d[i] = a[i] + I * b[i];FFT(d, N, op);d[N] = d[0];if(op == 1){for(int i = 0; i < N; i++){a[i] = (d[i] + d[N - i].conj()) / 2;b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;}} else {for(int i = 0; i < N; i++){a[i] = cp(d[i].x, 0);b[i] = cp(d[i].y, 0);}}d[N] = cp(0, 0);}void mult(int *a, int *b, int M){int N = 1, len = 0;while(N < M) N <<= 1, len++;for(int i = 0; i < N; i++)wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));for(int i = 0; i < N; i++){ka[i].x = a[i] >> 15;kb[i].x = b[i] >> 15;ta[i].x = a[i] & 32767;tb[i].x = b[i] & 32767;}D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1);for(int i = 0; i < N; i++){kk[i] = ka[i] * kb[i];kt[i] = ka[i] * tb[i] + ta[i] * kb[i];tt[i] = ta[i] * tb[i];}D_FFT(tt, kk, N, -1); FFT(kt, N, -1);for(int i = 0; i < N; i++){tt[i] = tt[i] / N;kt[i] = kt[i] / N;kk[i] = kk[i] / N;}}}PLE;void main() {scanf("%d%d%d", &n, &m, &p); n++; m++;for(int i = 0; i < n; i++) scanf("%d", &a[i]),a[i] = a[i] % p;for(int i = 0; i < m; i++) scanf("%d", &b[i]),b[i] = b[i] % p;PLE.mult(a, b, n + m - 1);for(int i = 0; i < n + m - 1; i++)printf("%lld ",(((((LL)round(kk[i].x)) % p) << 30) + ((((LL)round(kt[i].x)) % p) << 15) + ((LL)round(tt[i].x)) % p) % p);} }

模板：struct polynie {CP getw(int m, int k, int op) {return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m));}int wz[MAXN];CP A[MAXN], B[MAXN], C[MAXN];void FFT(CP *a, int N, int op) {rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);for(int l = 2; l <= N; l <<= 1) {int mid = l >> 1;CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid));for(int i = 0; i < N; i += l) {w = CP(1, 0);rop(j, 0, mid) {x = a[i + j];y = w * a[i + j + mid];a[i + j] = x + y;a[i + j + mid] = x - y;w = w * wn;}}}}void mult(CP *a, CP *b, CP *c, int M) {int N = 1, len = 0;while(N < M) N <<= 1, len++;rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));FFT(a, N, 1); FFT(b, N, 1);rop(i, 0, N) c[i] = a[i] * b[i];FFT(c, N, -1);rop(i, 0, N) c[i].x = c[i].x / N, c[i].y = c[i].y / N;}void blue_stein(CP *a, int M, int op) {int M2 = M << 1;memset(A, 0, sizeof(A));memset(B, 0, sizeof(B));memset(C, 0, sizeof(C));rop(i, 0, M) A[i] = a[i] * getw(M2, 1ll * i * i % M2, op);rop(i, 0, M2) B[i] = getw(M2, 1ll * (i - M) * (i - M) % M2, -op);mult(A, B, C, M2 + M - 1);rop(i, 0, M) a[i] = C[i + M] * getw(M2, 1ll * i * i % M2, op);if(op == -1) rop(i, 0, M) a[i].x = a[i].x / M, a[i].y = a[i].y / M;} }PLE;

多項式全集之三 多項式求逆:

多項式求逆：

1、問題描述：

已知

，且 ，求

2、推導過程：

設

由于

所以

那么

兩邊平方，得：

由于

的第 項為

一定有一項 ，所以

兩邊同乘

，得：

喜聞樂見的模板：

namespace INV{typedef long long LL;int n, a[300005], b[300005];const int mod = 998244353;const int g = 3189;int qpow(int a, int b){int ans = 1;while(b){if(b & 1) ans = 1ll * ans * a % mod;a = 1ll * a * a % mod;b >>= 1;}return ans;}struct p_l_e{int wz[300005], i_c[300005];void NTT(int *a, int N, int op){for(int i = 0; i < N; i++)if(i < wz[i]) swap(a[i], a[wz[i]]);for(int le = 2; le <= N; le <<= 1){int mid = le >> 1, wn = qpow(g, (mod - 1) / le);if(op == -1) wn = qpow(wn, mod - 2);for(int i = 0; i < N; i += le){LL w = 1; int x, y;for(int j = 0; j < mid; j++){x = a[i + j];y = w * a[i + j + mid] % mod;a[i + j] = (x + y) % mod;a[i + j + mid] = (x - y + mod) % mod;w = w * wn % mod;}}}}int init(int M){int N = 1, len = 0;while(N < M) N <<= 1, len++; for(int i = 0; i < N; i++)wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));return N;}void INV(int *a, int *b, int deg){if(deg == 1){b[0] = qpow(a[0], mod - 2); return;}INV(a, b, (deg + 1) >> 1);int N = init(deg + deg - 1);for(int i = 0; i < deg; i++) i_c[i] = a[i];for(int i = deg; i < N; i++) i_c[i] = 0;NTT(b, N, 1);NTT(i_c, N, 1);for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * (2 - 1ll * b[i] * i_c[i] % mod + mod) % mod;NTT(b, N, -1);int t = qpow(N, mod - 2);for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;for(int i = deg; i < N; i++) b[i] = 0;}}PLE;void main(){scanf("%d", &n);for(int i = 0; i < n; i++) scanf("%d", &a[i]);PLE.INV(a, b, n);for(int i = 0; i < n; i++) printf("%d ",b[i]);} }

多項式全集之四 多項式ln:

1、做法：

設

兩邊求導得

積分回去即可。

2、應用：

這個的組合意義是：無序組合。

設

， 表示一些東西，那么這些東西有序組合的方案數為

而無序組成的方案數為：

如果無序組合方案數好求，那么求

就能得到 。

例題

喜聞樂見的代碼：

多項式

:namespace PLE_ln{struct polyme {int li[SZ], wz[SZ];void NTT(int *a, int N, int op) {rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);for(int l = 2; l <= N; l <<= 1) {int mid = l >> 1;int x, y, w, wn = qpow(g, (mod - 1) / l);if(op) wn = qpow(wn, mod - 2);for(int i = 0; i < N; i += l) {w = 1;for(int j = 0; j < mid; ++j) {x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod;a[i + j] = (x + y) % mod;a[i + j + mid] = (x - y + mod) % mod;w = 1ll * w * wn % mod;}}}}void qd(int *a, int *b, int n) {rop(i, 0, n) b[i] = 1ll * a[i + 1] * (i + 1) % mod;}void jf(int *a, int *b, int n) {rop(i, 1, n) b[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;}void mult(int *a, int *b, int *c, int M) {int N = 1, len = 0;while(N < M) N <<= 1, len ++;rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));NTT(a, N, 0); NTT(b, N, 0);rop(i, 0, N) c[i] = 1ll * a[i] * b[i] % mod;NTT(c, N, 1);int t = qpow(N, mod - 2);rop(i, 0, N) c[i] = 1ll * c[i] * t % mod;}void inv(int *a, int *b, int deg) {if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;}inv(a, b, (deg + 1) >> 1);rop(i, 0, deg) li[i] = a[i];int N = 1, len = 0;while(N < deg + deg - 1) N <<= 1, len ++;rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));rop(i, deg, N) li[i] = 0;NTT(li, N, 0); NTT(b, N, 0);rop(i, 0, N) b[i] = 1ll * b[i] * (2 - 1ll * li[i] * b[i] % mod + mod) % mod;NTT(b, N, 1);int t = qpow(N, mod - 2);for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;rop(i, deg, N) b[i] = 0;}}PLE;int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n;void main() {scanf("%d", &n);rop(i, 0, n) scanf("%d", &a[i]);PLE.qd(a, da, n);PLE.inv(a, ia, n);PLE.mult(ia, da, dla, n + n - 1);PLE.jf(dla, la, n);rop(i, 0, n) printf("%d ", la[i]);} }

與作者交流：https://ac.nowcoder.com/discuss/175409

更多算法和題解：https://ac.nowcoder.com/acm/contest/discuss

總結

以上是生活随笔為你收集整理的蝴蝶优化算法_算法｜FFT基础及各种常数优化，5万字笔记：公式推导+代码模板...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。