学习笔记：FFT

复数

复数的概念

我们把形如 \(a+b\mathrm{i}(a, b\in \mathbf{R})\) 的数叫做复数，其中 \(\mathrm{i}\) 叫做虚数单位，\(\mathrm{i} = \sqrt {-1}\)。

全体复数构成的集合 \(\mathbf{C}=\{a+b\mathrm{i} \mid a,b\in \mathbf{R}\}\) 叫做复数集。

复数通常用字母 \(z\) 表示，即 \(z=a+b\mathrm{i}(a, b\in \mathbf{R})\)，其中 \(a\) 与 \(b\) 分别是复数 \(z\) 的实部与虚部，分别记作 \(\Re(z)\) 和 \(\Im(z)\)。

可以建立复平面来表示复数，即在直角坐标系 \(xOy\) 中，可以用点 \(Z(a,b)\) 来表示复数 \(z=a+b\mathrm{i}\)，其中 \(x\) 轴叫做实轴，\(y\) 轴叫做虚轴。复数 \(z\) 也与向量 \(\vec{OZ}\) 一一对应。向量 \(\vec{OZ}\) 的模叫做复数 \(z\) 的模或绝对值，记作 \(|z|\) 或 \(|a+b\mathrm{i}|\)，即 \(|z|=|a+b\mathrm{i}|=\sqrt{a^2+b^2}\)。

一般为了方便起见，通常把复数 \(z=a+b\mathrm{i}\) 说成点 \(Z\) 或说成向量 \(\vec{OZ}\)。

一般地，当两个复数实部相等，虚部互为相反数时，这两个复数叫做互为共轭复数。复数 \(z\) 的共轭复数用 \(\overline{z}\) 表示，即如果 \(z=a+b\mathrm{i}\)，那么 \(\overline{z}=a-b\mathrm{i}\)。

复数的四则运算

复数的加法

设 \(z_1=a+b\mathrm{i}\)，\(z_2=c+d\mathrm{i}\)，则 \(z_1+z_2=(a+c)+(b+d)\mathrm{i}\)。

满足交换律和结合律。

复数的减法

设 \(z_1=a+b\mathrm{i}\)，\(z_2=c+d\mathrm{i}\)，则 \(z_1-z_2=(a-c)+(b-d)\mathrm{i}\)。

复数的乘法

设 \(z_1=a+b\mathrm{i}\)，\(z_2=c+d\mathrm{i}\)，则 \(z_1z_2=(ac-bd)+(ad+bc)\mathrm{i}\)。

满足交换律，结合律和分配律。

复数的除法

设 \(z_1=a+b\mathrm{i}\)，\(z_2=c+d\mathrm{i}\)，则 \(z_1\div z_2=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}\mathrm{i}\)。将 \(\frac{z_1}{z_2}\) 上下同时乘 \(\overline{z_2}\) 即可。

复数的三角表示

一般地，任何一个复数 \(z=a+b\mathrm{i}\) 都可以表示为 \(r(\cos \theta + \mathrm{i}\sin\theta)\)，叫做复数 \(z\) 的三角表示式，其中 \(r\) 是 \(z\) 的模，\(\theta\) 是以实轴非负半轴为始边，射线 \(OZ\) 为终边的角，叫做 \(z\) 的辐角。

一个复数对应着无限多个辐角。我们规定在 \(0\le\theta<2\pi\)（有些地方好像定义为 \((-\pi,\pi]\)）的范围内的 \(z\) 的辐角 \(\theta\) 的值为辐角的主值，通常记作 \(\arg z\)。

\(z_1z_2 = r_1(\cos\theta_1+\mathrm{i}\sin\theta_1)\cdot r_2(\cos\theta_2+\mathrm{i}\sin\theta_2) = r_1r_2[\cos(\theta_1+\theta_2)+\mathrm{i}\sin(\theta_1+\theta_2)]\)，也就是两个复数相乘，模相乘，辐角相加。

\(\frac{z_1}{z_2}=\frac{r_1(\cos\theta_1+\mathrm{i}\sin\theta_1)}{r_2(\cos\theta_2+\mathrm{i}\sin\theta_2)} = \frac{r_1}{r_2}[\cos(\theta_1-\theta_2)+\mathrm{i}\sin(\theta_1-\theta_2)]\)，即两个复数相除，模相除，辐角相减。

欧拉公式

对于任意实数 \(x\) 有 \(\mathrm{e}^{\mathrm{i}x}=\cos x + \mathrm{i}\sin x\)。

复指数函数

对于复数 \(z=a+b\mathrm{i}\)，定义复指数函数为 \(\exp z = \mathrm{e}^z = \mathrm{e}^a(\cos b + \mathrm{i}\sin b)\)。

• \(|\exp z| = \exp a > 0\)。
• \(\arg \exp z = b\)。
• \(\exp(z_1+z_2)=\exp(z_1)\exp(z_2)\)。
• \(\exp z\) 以 \(2\pi \mathrm{i}\) 为最小正周期。

复三角函数

对于复数 \(z=a+b\mathrm{i}\)，定义复三角函数 \(\cos z = \frac{\exp(\mathrm{i}z) + \exp(-\mathrm{i}z)}{2}\)，\(\sin z = \frac{\exp(\mathrm{i}z) + \exp(-\mathrm{i}z)}{2\mathrm{i}}\)。

• 若取 \(z\in \mathbf{R}\)，则由欧拉公式，有，\(\cos z = \Re(\mathrm{e}^{\mathrm{i}z})\)，\(\sin z = \Im(\mathrm{e}^{\mathrm{i}z})\)。
• 复三角函数的正余弦函数均以 \(2\pi\) 为最小正周期。

复数的表示形式

代数形式：\(z=a+b\mathrm{i}\)。

三角形式：\(z=r(\cos\theta + \mathrm{i}\sin\theta)\)。

指数形式：\(z=r\exp(i\theta)\)。

单位根

若复数 \(z\) 满足 \(z^n=1\)，则称 \(z\) 为 \(n\) 次单位根。

设 \(\omega_n=\exp \frac{2\pi\mathrm{i}}{n}=\cos(\frac{2\pi}{n}) + \mathrm{i}\sin(\frac{2\pi}{n})\)（即辐角为 \(\frac{2\pi}{n}\) 的单位复数），则所有 \(n\) 次单位根 \(z\) 的集合可以表示为 \(\{\omega_n^{k}\mid k=0,1,\cdots,n-1\}\)，即 \(z_k = \omega_{n}^{k}\)。

每个单位根 \(z_k\) 对应复平面单位圆上一个点，其辐角为 \(\frac{2k\pi}{n}\)。

所有单位根将单位圆分成 \(n\) 等份，构成正 \(n\) 边形的顶点。

性质：

• \(\omega_{n}^{n}=1\)。
• \(\omega_{n}^{k}=\exp \frac{2k\pi\mathrm{i}}{n} = \cos(\frac{2k\pi}{n})+\mathrm{i}\sin(\frac{2k\pi}{n}) = \cos(\frac{2kx\pi}{nx})+\mathrm{i}\sin(\frac{2kx\pi}{nx}) = \omega_{nx}^{kx}\)。
• \(\omega_{2n}^{k+n}=-\omega_{2n}^{k}\)。

本原单位根

称集合 \(\{\omega_{n}^{k}\mid 0\le k<n,\gcd(n,k)=1\}\) 中的元素为本原单位根。

对于任何一个本原单位根 \(\omega_n\)，它的 \(k(0\le k<n)\) 次幂互不相同。从互质和整除的角度出发容易证明。

快速傅里叶变换 FFT

一种能在 \(O(n\log n)\) 的时间内计算两个给定的 \(n\) 次多项式的乘法的算法。

离散傅里叶变换 DFT & 离散傅里叶逆变换 IDFT

求出一个 \(n\) 次多项式在 \(n\) 个 \(n\) 次单位根下的点值的过程称为离散傅里叶变换 DFT。

而将这些点值重新插值成系数表示法的过程称为离散傅里叶逆变换 IDFT。

DFT：\(F_k = \sum\limits_{i=0}^{n-1}G_i\omega_{n}^{ik}\)。

IDFT：\(G_k=\frac{1}{n} \sum\limits_{i=0}^{n-1}F_i\omega_{n}^{-ik}\)。（下面有证明）

若要用 FFT 计算多项式乘法 \(C(x)=A(x)\cdot B(x)\)。其基本思想就是，在 \(O(n\log n)\) 时间内进行两遍 DFT 分别求出 \(A(x)\) 和 \(B(x)\) 的点值，再 \(O(n)\) 将点值相乘合并为 \(C(x)\) 的点值，最后在 \(O(n\log n)\) 时间内进行 IDFT 插值成 \(C(x)\) 的系数表示法。

快速傅里叶变换 FFT

FFT 算法可以 \(O(n\log n)\) 实现 DFT。

设 \(n-1\) 次多项式 \(A(x)=a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1}\)，

为了方便处理，以下我们将 \(A(x)\) 项数保证为 \(2\) 的整数次幂，若次数不够，则补系数 \(0\) 即可。即 \(A(x)\) 是不超过 \(n-1\) 次多项式，且 \(n\) 是 \(2\) 的整数次幂。

考虑如何快速计算其在 \(n\) 个 \(n\) 次单位根下的点值，即 \(A(\omega_{n}^0),A(\omega_{n}^1),\dots,A(\omega_{n}^{n-1})\)，

设，

\[A_1(x) = a_0 + a_2x + a_4x^2 + \cdots + a_{n-2}x^{\frac{n}{2}-1}， \\ A_2(x) = a_1 + a_3x + a_5x^2 + \cdots + a_{n-1}x^{\frac{n}{2}-1}， \]

则 \(A(x) = A_1(x^2) + xA_2(x^2)\)，

设 \(0\le k\le \frac{n}{2}-1\)，则，

\[\begin{aligned} A(\omega_{n}^{k}) &= A_1(\omega_{n}^{2k}) + \omega_{n}^{k}A_2(\omega_{n}^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) + \omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k}) \end{aligned} \]

又有，

\[\begin{aligned} A(\omega_{n}^{k+\frac{n}{2}}) &= A_1(\omega_{n}^{2k+n}) + \omega_{n}^{k+\frac{n}{2}}A_2(\omega_{n}^{2k+n}) \\ &= A_1(\omega_{n}^{2k}) - \omega_{n}^{k}A_2(\omega_{n}^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) - \omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k}) \end{aligned} \]

那么我们只需要求出 \(A_1(\omega_{\frac{n}{2}}^{k})\) 和 \(A_2(\omega_{\frac{n}{2}}^{k})\)，就可以知道 \(A(\omega_{n}^{k})\) 和 \(A(\omega_{n}^{k+\frac{n}{2}})\)，

由于 \(k\) 和 \(k+\frac{n}{2}\) 取遍了 \([0,n-1]\)，于是我们可以利用 \(A_1(x)\) 和 \(A_2(x)\) 求出 \(A(x)\) 在 \(n\) 个 \(n\) 次单位根下的点值，

由于 \(A_1(x)\) 和 \(A_2(x)\) 的项数都是 \(\frac{n}{2}\)，即为 \(A(x)\) 问题规模的一半的子问题，于是可以在 \(T(n)=2T(\frac{n}{2})+O(n)=O(n\log n)\) 的复杂度内求出 \(A(x)\) 的 \(n\) 个点值。

快速傅里叶逆变换 IFFT

IFFT 可以用 FFT 的方法实现，可以 \(O(n\log n)\) 实现 IDFT。

我们已知 \(y_i=A(\omega_{n}^{i})\)，现在要快速求出 \(A(x)\) 的系数 \(a_0,a_1,\dots,a_{n-1}\)，

考虑构造如下多项式，

\[F(x) = \sum\limits_{i=0}^{n-1} y_ix^i \]

设 \(b_i=\omega_{n}^{-i}\)，则，

\[\begin{aligned} F(b_k) &= \sum\limits_{i=0}^{n-1}A(\omega_{n}^{i})\omega_{n}^{-ki} \\ &= \sum\limits_{i=0}^{n-1}\omega_{n}^{-ki}\sum\limits_{j=0}^{n-1}a_{j}\omega_{n}^{ij} \\ &= \sum\limits_{j=0}^{n-1}a_j\sum\limits_{i=0}^{n-1}\omega_{n}^{(j-k)i} \end{aligned} \]

记 \(S(\omega_{n}^{x})=\sum\limits_{i=0}^{n-1}\omega_{n}^{xi}\)，

当 \(x\bmod n = 0\) 时，\(S(\omega_{n}^{x})=n\)。

当 \(x\bmod n \neq 0\) 时，

\[\begin{aligned} S(\omega_{n}^{x}) &= \sum\limits_{i=0}^{n-1}\omega_{n}^{xi} & (1) \\ \omega_{n}^{x}S(\omega_{n}^{x}) &= \sum\limits_{i=1}^{n}\omega_{n}^{xi} & (2) \\ \end{aligned} \]

\((2)-(1)\) 得，

\[S(\omega_{n}^{x}) = \frac{\omega_{n}^{xn}-\omega_{n}^{0}}{\omega_{n}^{x}-1} = 0 \]

因此，

\[\begin{aligned} F(b_k)&=a_k\cdot n \\ a_k &= \frac{F(\omega_{n}^{-k})}{n}=\frac{F(\omega_{n}^{n-k})}{n} \end{aligned} \]

同时也证明了 IDFT。

于是我们相当于要对构造多项式 \(F(x)\) 求出 \(n\) 次单位根下的点值 \(F(\omega_{n}^{k})\)，然后除以 \(n\) 得到 \(A(x)\) 的系数 \(a_{n-k}\)。

或者可以 FFT 直接求出 \(F(\omega_{n}^{-k})\) 的点值，再除以 \(n\)。由于 \(\omega_{n}^{-k}=(\omega_{n}^{-1})^k\)，相当于 FFT 证明过程中多带了一个 \(-1\) 次幂，不影响 FFT 的推导。

这里可以用数值分析的方法证明除以 \(n\) 四舍五入的答案不会受精度影响。

解释

FFT 为什么要选择单位根求值？

1. 首先，初始的单位根必须选择本原单位根，保证了其 \(0\) 到 \(n-1\) 次幂一定互不相同，因此就能确定 \(n-1\) 次多项式。这也说明了 IFFT 时选择 \(\omega_n^{-1}=\omega_n^{n-1}\) 的正确性。

2. 其次，若我们随便选择一个整数 \(v\)，将 \(v^0,v^1,v^2,\cdots,v^{n-1}\) 代入求值，那么一样也可以推出 \(A(v^k)=A_1(v^{2k})+v^kA_2(v^{2k})\)，但是只有在 \(k\in[0,\frac{n}{2}-1]\) 时成立。

   单位根的好处就是其利用了以下几个性质：

   1. \(\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}\)。
   2. \(\omega_{n}^{k+n}=\omega_{n}^{k}\)。
   3. \(\omega_{mn}^{mk}=\omega_{n}^{k}\)。

NTT 即选择了与单位根有类似性质的原根来代入求值，相对于 FFT 的优点就是常数较小，且没有精度误差。但是 NTT 只在对一些特殊质数取模才有效。

代码实现

模板题：LGP3803 【模板】多项式乘法（FFT）

分治法

用一个函数可以同时实现 FFT 和 IFFT。

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
using ULL = unsigned long long;
using LD = long double;
using PII = pair<int, int>;

const int N = 1e6 + 5;
const double PI = acos(-1.0);

int n, m;

void FFT(vector<complex<double>>& a, int n, int o) {
    if (n == 1) return;
    vector<complex<double>> a1(n / 2), a2(n / 2);
    for (int i = 0; i < n; i += 2) a1[i / 2] = a[i], a2[i / 2] = a[i + 1];
    FFT(a1, n / 2, o), FFT(a2, n / 2, o);
    complex<double> p(1, 0), w(cos(2 * PI / n), o * sin(2 * PI / n));
    for (int i = 0; i < n / 2; i++) {
        a[i] = a1[i] + p * a2[i];
        a[i + n / 2] = a1[i] - p * a2[i];
        p *= w;
    }
}

int main() {
    cin >> n >> m;
    int tot = 1;
    while (tot < n + m + 1) tot <<= 1;
    vector<complex<double>> a(tot), b(tot);
    for (int i = 0; i <= n; i++) {
        double x;
        cin >> x;
        a[i].real(x);
    }
    for (int i = 0; i <= m; i++) {
        double x;
        cin >> x;
        b[i].real(x);
    }
    FFT(a, tot, 1);
    FFT(b, tot, 1);
    vector<complex<double>> c(tot);
    for (int i = 0; i < tot; i++) c[i] = a[i] * b[i];
    FFT(c, tot, -1);
    for (int i = 0; i <= n + m; i++) cout << (int)round(c[i].real() / tot) << ' ';
    cout << '\n';
    return 0;
}

蝶形变换法

分治法实现的常数非常大，且空间复杂度是 \(O(n\log n)\) 的，于是有了更高效的蝶形变换法。

考虑 FFT 的时候，每次分治都要把奇偶位的系数分开，那么会不会存在什么规律？

假设 \(n-1\) 次 \(n\) 项多项式 \(A(x)\) 的系数有 \(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\)，\(n=8\)，在 FFT 分治到最后的时候，系数位置会变成 \(a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7\)，

我们发现 \(a_i\) 位置上的系数，最后变成 \(a_j\)，满足 \(j\) 是 \(i\) 在补全 \(m=\log_2n\) 位二进制后的二进制翻转。例如 \((000)_2\to(000)_2\)，\((001)_2\to(100)_2\)，\((110)_2\to (011)_2\)。

证明：

考虑系数 \(a_x\)。下证分治过后，\(a_x\) 的位置会变成 \(y\)，\(y\) 是 \(x\) 在补全 \(m\) 位二进制后的二进制翻转。

设 \(x\) 的二进制第 \(i\) 位是 \(x_i\)，\(i\) 从 \(0\) 开始。

分治一共会进行 \(m\) 层，第 \(m\) 层只剩下一个系数，也就是区间长度为 \(1\)。

从 \(1\) 到 \(m-1\) 考虑每一层，我们模拟一遍递归的过程。

若 \(x\) 最低位是 \(1\)，则下一层 \(x\to \frac{n}{2}+\frac{x-1}{2}\)，但是下一层之后我们为了方便计算 \(x\) 的目标位置，因此不能把 \(\frac{n}{2}\) 带下去，于是记一个 \(k\)（第一层之前 \(k=0\)），然后 \(k\to k+\frac{n}{2}\)，\(x\to \frac{x-1}{2}\)，\(n\to \frac{n}{2}\)。

若 \(x\) 最低位是 \(0\)，则下一层 \(k\to k+0\)，\(x\to \frac{x}{2}\)， \(n\to \frac{n}{2}\)。

那么每一层相当于把 \(x\) 当前最低位，以 \(x\) 二进制位中线为对称轴翻转，计入变量 \(k\)，然后把 \(x\) 最低位去掉（除以 \(2\)）。

那么这样递归过后，最终 \(x\) 的值变为 \(0\)，然后要加上 \(k\) 才是 \(x\) 在原序列中的位置，容易看出 \(k\) 就是一开始的 \(x\) 的二进制翻转。

证毕。

于是我们可以提前交换好系数 \(a_i\) 的位置，然后自底向上合并，合并的时候可以在原数组上进行覆写，那么空间复杂度就降低到了 \(O(n)\)，且常数比分治法的递归小了很多。

如何求一个数的二进制翻转呢？这里可以 \(O(n)\) 递推实现。

设 \(\mathrm{rev}(x)\) 表示 \(x\) 的二进制翻转，显然 \(\mathrm{rev}(0)=0\)，

假设我们要求 \(\mathrm{rev}(x)\)，且 \(\mathrm{rev}(\lfloor\frac{x}{2}\rfloor)\) 已经求得，于是有，

\[\mathrm{rev}(x)=\lfloor\frac{\mathrm{rev}(\lfloor\frac{x}{2}\rfloor)}{2}\rfloor + [x\bmod 2=1]\cdot \frac{n}{2} \]

例如 \(n=32\)，\((01101)_2\to (00110)_2 \to (01100)_2\to (00110)_2 \to(10110)\)。

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
using ULL = unsigned long long;
using LD = long double;
using PII = pair<int, int>;

const int N = 4e6 + 5;
const double PI = acos(-1.0);

int n, m;
complex<double> a[N], b[N], c[N];
int rev[N];

void FFT(complex<double> *a, int n, int o) {
    rev[0] = 0;
    for (int i = 1; i < n; i++) rev[i] = (rev[i >> 1] >> 1) + (i & 1) * n / 2;
    for (int i = 0; i < n; i++) if (rev[i] < i) swap(a[i], a[rev[i]]);
    for (int len = 2; len <= n; len <<= 1) {
        for (int l = 0; l + len - 1 < n; l += len) {
            int mid = len >> 1;
            complex<double> w(cos(2 * PI / len), o * sin(2 * PI / len)), p(1, 0);
            for (int i = 0; i < mid; i++) {
                auto x = a[l + i], y = a[l + i + mid];
                a[l + i] = x + p * y;
                a[l + i + mid] = x - p * y;
                p *= w;
            }
        }
    }
    if (o == -1) for (int i = 0; i < n; i++) a[i] /= n;
}

int main() {
    ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
    cin >> n >> m;
    double x;
    for (int i = 0; i <= n; i++) {
        cin >> x;
        a[i].real(x);
    }
    for (int i = 0; i <= m; i++) {
        cin >> x;
        b[i].real(x);
    }
    int tot = 1;
    while (tot < n + m + 1) tot <<= 1;
    FFT(a, tot, 1);
    FFT(b, tot, 1);
    for (int i = 0; i < tot; i++) c[i] = a[i] * b[i];
    FFT(c, tot, -1);
    for (int i = 0; i <= n + m; i++) cout << (int)round(c[i].real()) << ' ';
    cout << '\n';
    return 0;
}

三步变两步

设多项式 \(F(x) = A(x) + B(x)\mathrm{i}\)，则 \(F(x)^2 = A(x)^2 - B(x)^2 + 2\cdot A(x)\cdot B(x) \cdot \mathrm{i}\)，于是 \(A(x)\cdot B(x) = \frac{\Im(F(x))}{2}\)。

于是若要求 \(A(x)\cdot B(x)\) 的系数，可以直接求出 \(F(x)^2\)，再取 \(F(x)^2\) 每一项系数的虚部除以 \(2\) 即可，只需两步 FFT。

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
using ULL = unsigned long long;
using LD = long double;
using PII = pair<int, int>;

const int N = 4e6 + 5;
const double PI = acos(-1.0);

int n, m;
complex<double> a[N];
int rev[N];

void FFT(complex<double> *a, int n, int o) {
    rev[0] = 0;
    for (int i = 1; i < n; i++) rev[i] = (rev[i >> 1] >> 1) + (i & 1) * n / 2;
    for (int i = 0; i < n; i++) if (rev[i] < i) swap(a[i], a[rev[i]]);
    for (int len = 2; len <= n; len <<= 1) {
        for (int l = 0; l + len - 1 < n; l += len) {
            int mid = len >> 1;
            complex<double> w(cos(2 * PI / len), o * sin(2 * PI / len)), p(1, 0);
            for (int i = 0; i < mid; i++) {
                auto x = a[l + i], y = a[l + i + mid];
                a[l + i] = x + p * y;
                a[l + i + mid] = x - p * y;
                p *= w;
            }
        }
    }
    if (o == -1) for (int i = 0; i < n; i++) a[i] /= n;
}

int main() {
    ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
    cin >> n >> m;
    double x;
    for (int i = 0; i <= n; i++) {
        cin >> x;
        a[i].real(x);
    }
    for (int i = 0; i <= m; i++) {
        cin >> x;
        a[i].imag(x);
    }
    int tot = 1;
    while (tot < n + m + 1) tot <<= 1;
    FFT(a, tot, 1);
    for (int i = 0; i < tot; i++) a[i] *= a[i];
    FFT(a, tot, -1);
    for (int i = 0; i <= n + m; i++) cout << (int)round(a[i].imag() / 2) << ' ';
    cout << '\n';
    return 0;
}

各优化时间对比

题目

P1919 【模板】高精度乘法 | A*B Problem 升级版

相当于给了两个多项式在 \(x=10\) 位置的点值，求多项式相乘后在 \(x=10\) 位置的点值。

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
using ULL = unsigned long long;
using LD = long double;
using PII = pair<int, int>;

const int N = 4e6 + 5;
const double PI = acos(-1.0);

string s, t;
complex<double> a[N];
int rev[N], ans[N];

void FFT(complex<double> *a, int n, int o) {
    rev[0] = 0;
    for (int i = 1; i < n; i++) rev[i] = (rev[i >> 1] >> 1) + (i & 1) * n / 2;
    for (int i = 0; i < n; i++) if (rev[i] < i) swap(a[i], a[rev[i]]);
    for (int len = 2; len <= n; len <<= 1) {
        for (int l = 0; l + len - 1 < n; l += len) {
            int mid = len >> 1;
            complex<double> w(cos(2 * PI / len), o * sin(2 * PI / len)), p(1, 0);
            for (int i = 0; i < mid; i++) {
                auto x = a[l + i], y = a[l + i + mid];
                a[l + i] = x + p * y;
                a[l + i + mid] = x - p * y;
                p *= w;
            }
        }
    }
}

int main() {
    cin >> s >> t;
    int n = s.size(), m = t.size();
    for (int i = n - 1; i >= 0; i--) a[i].real(s[n - i - 1] - '0');
    for (int i = m - 1; i >= 0; i--) a[i].imag(t[m - i - 1] - '0');
    int tot = 1;
    while (tot < n + m - 1) tot <<= 1;
    FFT(a, tot, 1);
    for (int i = 0; i < tot; i++) a[i] *= a[i];
    FFT(a, tot, -1);
    int len = n + m - 1;
    for (int i = 0; i < len; i++) {
        ans[i] += round(a[i].imag() / tot / 2);
        if (ans[i] >= 10) {
            ans[i + 1] += ans[i] / 10;
            ans[i] %= 10;
            len += (i == len - 1);
        }
    }
    while (len > 1 && ans[len - 1] == 0) len--;
    for (int i = len - 1; i >= 0; i--) cout << ans[i];
    cout << '\n';
    return 0;
}

HDU4609 3-idiots

设 \(f(i)\) 表示相加起来 \(=i\) 的两根木棍的对数。

设 \(g(i)\) 表示长度 \(=i\) 的木棍的个数，即对 \(a_i\) 开一个桶。

把 \(g\) 看成多项式，\(g(i)\) 是 \(i\) 次项的系数，设 \(h=g^2\)，则，

\[h(i) = \sum\limits_{j=1}^{i} g(j)\cdot g(i-j) \]

\(h\) 用 FFT 计算即可。

那么 \(f(i) = h(i) - \sum\limits_{j=1}^{n} [a_j+a_j=i]\)，即减去两个重复选同一根的方案。

于是我们对 \(a_i\) 升序排序，枚举 \(i\in [1,n]\)，计算以 \(a_i\) 为三根木棍中的最长木棍，且三根木棍能构成三角形的方案数。

• 首先答案加上 \(\sum\limits_{j=a_i+1} f(j)\)，表示选择的另外两根木棍长度和 \(>a_i\)。
• 由于我们要保证两根最小的长度加起来大于第三根，所以剩下两根必须选择 \(\le a_i\) 的，并且不能重复选 \(a_i\)，因此要减去一些不合法的方案。
• 减去 \((n-i)\cdot (i-1)\)，表示一根选择了 \(\ge a_i\) 的，一根选择 \(\le a_i\) 的，但都没有重复选 \(a_i\) 的方案数。
• 减去 \(n-1\)，表示一根选择了 \(a_i\) 这根（注意是 \(a_i\) 这根，而不是长度 \(=a_i\)，也就是选重复了），另一根选了其他的方案数。
• 减去 \((n-i)\cdot (n-i-1) \cdot\frac{1}{2}\)，表示两个都选择了 \(\ge a_i\)，都没有重复选 \(a_i\) 的方案数。
• 两个都重复选 \(a_i\) 的方案数已经在前面计算 \(f\) 的时候减过了。

最后合法方案除以总方案就是概率了。

#include <bits/stdc++.h>
using namespace std;
using LL = long long;
using ULL = unsigned long long;
using LD = long double;
using PII = pair<int, int>;

const int N = 4e5 + 5, V = 2e5;
const double PI = acos(-1.0);

int n, a[N], tot;
LL b[N], f[N], sum[N];
complex<double> B[N];

namespace Poly {
int rev[N];

void FFT(complex<double> *a, int n, int o) {
    rev[0] = 0;
    for (int i = 1; i < n; i++) rev[i] = (rev[i >> 1] >> 1) + (i & 1) * n / 2;
    for (int i = 0; i < n; i++) if (rev[i] < i) swap(a[i], a[rev[i]]);
    for (int len = 2; len <= n; len <<= 1) {
        for (int l = 0; l + len - 1 < n; l += len) {
            int mid = len >> 1;
            complex<double> p(1, 0), w(cos(2 * PI / len), o * sin(2 * PI / len));
            for (int i = 0; i < mid; i++) {
                auto x = a[l + i], y = a[l + i + mid];
                a[l + i] = x + p * y;
                a[l + i + mid] = x - p * y;
                p *= w;
            }
        }
    }
}
} using Poly::FFT;

void Solve() {
    cin >> n;
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
        b[a[i]]++;
    }

    sort(a + 1, a + 1 + n);

    int m = a[n], tot = 1;
    while (tot < m + m + 1) tot <<= 1;

    // 计算 f = b * b
    for (int i = 0; i <= m; i++) B[i].real(b[i]);
    FFT(B, tot, 1);
    for (int i = 0; i < tot; i++) B[i] *= B[i];
    FFT(B, tot, -1);
    for (int i = 0; i < tot; i++) B[i] /= tot;
    for (int i = 0; i < tot; i++) f[i] = round(B[i].real());

    for (int i = 1; i <= n; i++) f[a[i] + a[i]]--; // 减去重复选的 a[i]+a[i]
    for (int i = 1; i < tot; i++) sum[i] = sum[i - 1] + f[i] / 2; // 除以2 是因为一对 a[i]+a[j] 会算两次

    LL ans = 0;
    // 枚举选取的最大长度是 a[i]
    for (int i = 1; i <= n; i++) {
        ans += sum[tot - 1] - sum[a[i]]; // 先计算上 加起来 >a[i] 的方案数
        // 因为 a[i] 的最大长度，因此我们要求剩余选的两个必须都 <=a[i]，且加起来 >a[i]
        ans -= (LL)(n - i) * (i - 1); // 一个选 >=a[i] 的，另一个选 <=a[i] 的，不合法
        ans -= n - 1; // 一个选 a[i]，另一个选不是 a[i] 的，不合法
        ans -= (LL)(n - i) * (n - i - 1) / 2; // 两个都选 >=a[i] 的，不合法
        // 两个都选 a[i] 已经在前面减过了
        // 剩下的就是两个都 <=a[i]，且加起来 >a[i] 的方案数
    }

    printf("%.7lf\n", (double)ans / ((LL)n * (n - 1) * (n - 2) / 6));
    
    // clear
    for (int i = 0; i < tot; i++) {
        b[i] = f[i] = sum[i] = 0;
        B[i] = {0, 0};
    }
}

int main() {
    int T;
    cin >> T;
    while (T--) Solve();
    return 0;
}