蝴蝶优化算法_算法｜FFT基础及各种常数优化，5万字笔记：公式推导+代码模板...

作者：中二攻子

链接：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