多项式

快速傅里叶变换（FFT）

前置知识：复数、单位根。

FFT 用于计算两个多项式的乘法。其时间复杂度为 $O(n\log n)$。

系数表示法和点值表示法

我们知道，$n+1$ 个不同的点 $(x,y)$ 可以确定一个唯一的 $n$ 次多项式。

那么，给出 $n+1$ 个不同的坐标来表示一个多项式，称为点值表示法；将多项式写为 $f(x)=\sum_{i=0}^n a_ix^i$，称为系数表示法

算法思想

我们有以下想法：把两个函数转换为点值表示法，将对应的点值纵坐标相乘，最后再将它转换为系数表示法，就得到了它们的乘积。

例如，我们可以求出 $x\in[1,n+m+1]$ 时 $F(x)$ 和 $G(x)$ 的值，将对应的点值纵坐标相乘，就得到了它们乘积的点值表示法。然后我们就可以解出答案的各项系数。

但是，例子的时间复杂度最优只能做到 $O(n^2)$。我们需要找出一些特殊的点值，使得上述过程的时间复杂度降低。这里我们使用单位根。

求出一个 $n$ 次多项式在每个 $n$ 次单位根下的点值，被称为离散傅里叶变换（Discrete Fourier transform，DFT）。相反，用这些点值求出多项式各项系数的过程，被称为离散傅里叶逆变换（Inverse Discrete Fourier transform，IDFT）。

DFT

假设我们有 $n$ 次多项式 $F(x)$，则我们要求出一个长度为 $n$ 的序列 $b$，其中 $b_i$ 表示 $F(x)$ 在 $n$ 次单位根的 $i$ 次幂下的点值，即：

$$b_i=\sum_{j=0}^{n-1} a_j\times (\omega_n^i)^j$$

暴力做是 $O(n^2)$ 的，我们考虑用快速傅里叶变换（Fast Fourier Transform，FFT）来优化这个过程。

FFT

我们考虑分治方法。

首先向 $A(x)$ 的最高次项添加 $0$，直到次数为 $2^k$。这样我们在分治的过程中可以将序列刚好分成两半。

设 $A$ 的各项系数为 $\{a_0,a_1,a_2,...,a_{n-1}\}$。

我们将所有偶数次项系数拿出来作为新函数 $A_0$，将所有奇数次项拿出来作为新函数 $A_1$，即：

$$A_0(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-2}x^{\frac{n}{2}-1}\\ A_1(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-1}x^{\frac{n}{2}-1}$$

容易得到 $A(x)=A_0(x^2)+xA_1(x^2)$。

我们先递归求出 $A_0$、$A_1$ 的变换结果，考虑如何由此推出 $A$ 的变换结果。

$$A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^{2k})$$

单位根有性质 $\omega_n^k=\omega_{2n}^{2k}$，可以推出 $\omega_n^{2k}=\omega_{\frac{n}{2}}^{k}$。带入得：

$$A(\omega_n^k)=A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})$$

此时，如果 $k<\frac{n}{2}$，则 $A_0(\omega_{\frac{n}{2}}^{k})$ 是 $A_0$ 变换结果的第 $k$ 项，$A_1(\omega_{\frac{n}{2}}^{k})$ 是 $A_1$ 变换结果的第 $k$ 项。由此我们可以求得 $A$ 变换结果的前 $\frac{n}{2}$ 项。

接下来我们求变换结果的后 $\frac{n}{2}$ 项。单位根有性质 $\omega_{n}^{k+\frac{n}{2}}=-\omega_n^k$。由 $A(x)=A_0(x^2)+xA_1(x^2)$ 可得，$A(-x)=A_0(x^2)-xA_1(x^2)$。所以可得：

$$A(\omega_{n}^{k+\frac{n}{2}})=A(-\omega_n^k)=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})$$

这样我们就求出了 $A$ 的变换结果。时间复杂度 $O(n \log n)$。

void FFT(int n, complex<double>* a)
{    
    if(n <= 1) return; 
    complex<double> a0[n / 2], a1[n / 2];
    for(int i = 0; i < n / 2; i++) a0[i] = a[i * 2], a1[i] = a[i * 2 + 1];
    FFT(n / 2, a0), FFT(n / 2, a1);
    for(int i = 0; i < n / 2; i++)
    {
        complex<double> x(cos(pi * 2 * i / n), sin(pi * 2 * i / n)); // 这个值是 w_n^k
        a[i] = a0[i] + x * a1[i];
        a[i + n / 2] = a0[i] - x * a1[i];
    }
}

逆变换

我们已经把两个多项式的乘积用点值表示法表示出来了，现在我们要把它转换为系数表示法。

假设我们已知离散傅里叶变换后的点值 ${b_i}$，我们想求多项式的系数 ${a_i}$。

先写公式：

$$a_i=\frac{1}{n}\sum_{j=0}^{n-1}b_j\times\omega_n^{-ij}$$

证明

由离散傅里叶变换的定义可得：

$$b_i=\sum_{j=0}^{n-1} a_j\times \omega_n^{ij}$$

将 $b_i$ 带入到 $\sum_{j=0}^{n-1}b_j\times\omega_n^{-ij}$（即上述公式去掉 $\frac{1}{n}$）：

$$\begin{aligned} \sum_{j=0}^{n-1}b_j\times\omega_n^{-ij} &=\sum_{j=0}^{n-1}\sum_{k=0}^{n-1} a_k\times \omega_n^{jk}\times \omega_n^{-ij}\\ &=\sum_{j=0}^{n-1}\sum_{k=0}^{n-1} a_k\times \omega_n^{j(k-i)}\\ &=\sum_{k=0}^n a_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)}\\ \end{aligned}$$

我们观察式子中的 $\sum_{j=0}^{n-1}\omega_n^{j(k-i)}$。此时分两种情况：

1、$i=k$

此时 $k-i=0$，所以：

$$\sum_{j=0}^{n-1}\omega_n^{j(k-i)}=\sum_{j=0}^{n-1}1=n$$

2、$i \neq k$

此时 $\sum_{j=0}^{n-1}\omega_n^{j(k-i)}$ 可以看成等比数列求和的形式：

$$\sum_{j=0}^{n-1}(\omega_n^{k-i})^j=\frac{1-(\omega_n^{k-i})^n}{1-\omega_n^{k-i}}$$

我么们观察分子中的 $(\omega_n^{k-i})^n$，它等于 $(\omega_n^n)^{k-i}$。单位根有性质 $\omega_n^n=1$，所以这个式子的值就是 $1$。所以：

$$\sum_{j=0}^{n-1}\omega_n^{j(k-i)}=0$$

综上，我们可以得到：

$$\sum_{j=0}^{n-1}b_j\times\omega_n^{-ij}=\sum_{k=0}^n a_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)}=\sum_{k=0}^n a_k[i=k]\times n=n\times a_i$$

所以：

$$\frac{1}{n}\sum_{j=0}^{n-1}b_j\times\omega_n^{-ij}=a_i$$

证毕。

递归实现

如果把逆变换式子中的 $\frac{1}{n}$ 去掉，发现它与变换的公式很相似。正变换求 $\omega_n^i$ 的点值，而逆变换相当于求 $\omega_n^{-i}$ 的点值。逆变换的推导过程和正变换十分相似，这里就不再展开。它们最终的实现过程也很相似，我们可以把变换和逆变换放在同一个函数中。

逆变换后要记住将所有项除以 $n$。

void FFT(int n, complex<double>* a, int p) // p = 1 为正变换，p = -1 为逆变换
{
    if(n <= 1) return;
    complex<double> a0[n / 2], a1[n / 2];
    for(int i = 0; i < n / 2; i++) a0[i] = a[i * 2], a1[i] = a[i * 2 + 1];
    FFT(n / 2, a0, p), FFT(n / 2, a1, p);
    complex<double> z(cos(pi * 2 / (p * n)), sin(pi * 2 / (p * n))), w(1, 0);
    for(int i = 0; i < n / 2; i++)
    {
        // 若 p = 1，则此时 w = w_n^i
        // 若 p = -1，则此时 w = w_n^{-i}
        // 不要在这里直接求 w 的值，因为 sin 和 cos 的计算速度较慢
        a[i] = a0[i] + w * a1[i];
        a[i + n / 2] = a0[i] - w * a1[i];
        w *= z;
    }
}

迭代实现

递归实现的时间复杂度虽然正确，但是递归、动态开数组使得它的常数较大。我们要找到更优的实现方法。

假设 $n=8$，我们来模拟它的递归过程：

第一层：0 1 2 3 4 5 6 7
第二层：0 2 4 6|1 3 5 7
第三层：0 4|2 6|1 5|3 7
第四层：0|4|2|6|1|5|3|7

递归过程的最底层是 0 4 2 6 1 5 3 7，我们来找它们的规律。我们先写出它们的二进制：

000 100 010 110 001 101 011 111

如果我们把这些二进制数翻转，就可以得到：

000 001 010 011 100 101 110 111

转换为十进制，它们就是：

0 1 2 3 4 5 6 7 8

现在我们发现，最底层第 $i$ 位上的数字，就是 $i$ 二进制下翻转后的数字。

所以我们要求一个数组 $rev_i$，表示 $i$ 二进制下翻转后的数字。

$i$ 的二进制，可以用 $\lfloor\frac{i}{2}\rfloor$ 二进制左移一位，再加上 $[i \mod 2=1]$ 计算出。翻转后，我们可以得到：$rev_i$ 可以用 $rev_{\lfloor\frac{i}{2}\rfloor}$ 右移一位，再加上 $[i \mod 2 = 1]\times(n>>1)$ 计算出。

int rev[N];
void make_rev(int n)
{
    for(int i = 1; i < n; i++)
    {
        rev[i] = rev[i >> 1] >> 1;
        if(i & 1) rev[i] |= (n >> 1);
    }
}

如果仔细思考，就会发现，最底层的变换，其实就是第 $i$ 位与第 $rev_i$ 位的数字交换位置。

如果知道了最底层的数字，我们逐层向上合并即可。

洛谷 P3803完整代码：

#include<bits/stdc++.h>
using namespace std;

const int N = 4000010;
int rev[N];
void make_rev(int n)
{
    for(int i = 1; i < n; i++)
    {
        rev[i] = rev[i >> 1] >> 1;
        if(i & 1) rev[i] |= (n >> 1);
    }
}
const double pi = acos(-1);
void FFT(int n, complex<double>* a, int p)
{
    make_rev(n);
    for(int i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]); // 将 a 按照最底层的数字排列
    for (int l = 2; l <= n; l <<= 1) //l：当前区间的长度
    {
        complex<double> z(cos(pi * 2 / (p * l)), sin(pi * 2 / (p * l)));
        for (int i = 0; i < n; i += l) //i：区间左端点
        {
            complex<double> w(1, 0);
            for (int j = i; j < i + l / 2; j++)
            {
                complex<double> x = a[j] + w * a[j + l / 2], y = a[j] - w * a[j + l / 2];
                a[j] = x, a[j + l / 2] = y;
                w *= z;
            }
        }
    }
}
complex<double> a[N], b[N], c[N];
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    int n, m, k = 1;
    cin >> n >> m;
    while(k < n + m + 1) k *= 2; // 找到第一个 >= n + m + 1 的二的整数次幂
    for (int i = 0; i <= n; i++)
    {
        int x;
        cin >> x;
        a[i] = complex<double>(x, 0);
    }
    for (int i = 0; i <= m; i++)
    {
        int x;
        cin >> x;
        b[i] = complex<double>(x, 0);
    }
    FFT(k, a, 1);
    FFT(k, b, 1);
    for (int i = 0; i < k; i++) c[i] = a[i] * b[i];
    FFT(k, c, -1);
    for (int i = 0; i <= n + m; i++) cout << (int)round(c[i].real() / k) << " ";
    return 0;
}

---

快速数论变换（NTT）

前置知识：FFT、原根。

在 FFT 中，我们用到了单位根处德点值，这就无法避免精度问题。在 NTT 中，我们放弃单位根，使用 原根。如果运算只涉及整数，或者结果要对 $998244353$ 这样的质数取模，就可以使用 NTT。

算法

我们可以取大质数 $p=998244353$，它的原根是 $g=3$。接下来的运算都在 $\mod p$ 的意义下。

对于 $n$（我们默认 $n$ 是二的整数次幂），我们设 $qn+1=p$，此时 $\omega_n=g^q$。也就是 $\omega_n=g^{\frac{p-1}{n}}$。

对于 $0\le i < n$，我们要求 $\omega_n^i$ 处的点值。

为什么原根可以代替单位根

我们思考为什么可以使用单位根：

1. 对于 $0\le i < n$，$\omega_n^i$ 互不相同。
2. $\omega_n^k=\omega_{2n}^{2k}$，用于变换。
3. $\omega_{n}^{k+\frac{n}{2}}=-\omega_n^k$，用于变换。
4. 对于 $k\neq 0$，$\sum_{i=0}^{n-1}\omega_n^{ik}=0$，用于逆变换。

如果此时 $\omega_n=g^q$，我们发现，上面的四条性质仍然成立！

性质 1：

原根有性质：若 $p$ 为质数，则 $g^i \mod p$，$0\le i < p-1$ 的结果两两不同。

所以 $1,g^{q},g^{2q},g^{3q}...$ 的结果互不相同。

性质 2：

$$\omega_n^k=g^{k\frac{p-1}{n}}=g^{2k\frac{p-1}{2n}}=\omega_{2n}^{2k}$$

性质 3：

由费马小定理得：

$$(\omega_n^{\frac{n}{2}})^2=\omega_n^n=g^{n\frac{p-1}{n}}=g^{p-1}\equiv 1\pmod p$$

所以 $\omega_n^{\frac{n}{2}}\equiv \pm 1\pmod p$。

根据性质一，$\omega_n^{\frac{n}{2}}\neq \omega_n^0=1$，所以 $\omega_n^{\frac{n}{2}}\equiv -1\pmod p$。

那么：

$$\omega_{n}^{k+\frac{n}{2}}=\omega_{n}^k\times \omega_n^{\frac{n}{2}}=-\omega_{n}^k$$

性质 4：

由等比数列求和公式得：

$$\sum_{i=0}^{n-1}\omega_n^{ik}=\frac{1-(\omega_n^k)^n}{1-\omega_n^k}$$

在性质 3 得推到中，我们推出了 $\omega_n^n=1$，所以分子中的 $(\omega_n^k)^n=(\omega_n^n)^k=1$，所以原式的值是 $0$。

实现

我们可以直接在 FFT 代码上修改，把 $w_n$ 从单位根换成原根就可以了。不要忘记在运算过程中随时取模。

const int p = 998244353;
void NTT(int n, int* a, int k)
{
    make_rev(n);
    for(int i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int l = 2; l <= n; l <<= 1)
    {
        int z = qpow(3, (p - 1) / l);
        if(k == -1) z = qpow(z, p - 2);
        for (int i = 0; i <= n - 1; i += l)
        {
            int w = 1;
            for (int j = i; j < i + l / 2; j++)
            {
                int x = (a[j] + 1ll * w * a[j + l / 2] % p) % p, y = (a[j] - 1ll * w * a[j + l / 2] % p + p) % p;
                a[j] = x, a[j + l / 2] = y;
                w = 1ll * w * z % p;
            }
        }
    }
}

其它常用的模数

摘自 OI-Wiki。

$$p=167772161=5\times2^{25}+1,g=3\\ p=469762049=7\times2^{26}+1,g=3\\ p=754974721=45\times2^{24}+1,g=11\\ p=998244353=119\times 2^{23}+1,g=3\\ p=1004535809=479\times 2^{21}+1,g=3$$

---

快速沃尔什变换（FWT）

FWT 用于解决下标位运算卷积问题，形式为：给定两个序列 $A_i$ 和 $B_i$，求出序列 $C$，满足：

$$C_i=\sum_{j\oplus k=i}A_jB_k$$

即 $A_jB_k$ 贡献在 $C$ 的第 $j\oplus k$ 位上。

其中 $\oplus$ 为位运算的某一种。通常为或、与、异或中的一种。

FWT 的核心思想是：对于一种运算，要找到一种构造方法。假设数组 $A$ 的构造结果是 $FWT[A]$，则需要满足 $FWT[C]_i=FWT[A]_i\times FWT[B]_i$。我们求出 $FWT[A]$ 和 $FWT[B]$，通过每一位的乘积求出 $FWT[C]$，最后通过逆变换求出 $C$。

或

构造

我们要求：

$$C_i=\sum_{j | k=i}A_jB_k$$

我们的构造方案是：

$$FWT[A]_i=\sum_{j|i=i}A_j$$

即 $i$ 所有子集对应位置之和。

构造正确性的证明

$$\begin{aligned} &~~~~~FWT[A]_i\times FWT[B]_i\\ &= \sum_{j|i=i}A_j\times\sum_{k|i=i}B_k\\ &= \sum_{j|i=i}~\sum_{k|i=i}A_jB_k\\ &= \sum_{(j|k)|i=i}A_jB_k\\ &= \sum_{x|i=i}~\sum_{j|k=x}A_jB_k\\ &= \sum_{x|i=i}C_x\\ &= FWT[C]_i \end{aligned}$$

变换（FWT）

按照式子直接做是 $O(n^2)$ 的，无法接受。

我们考虑分治。假设当前有长度为 $2^k$ 的序列 $A_{0\sim 2^k-1}$，它的前半段为 $A_0$，后半段是 $A_1$。假设我们已经递归求出了 $FWT[A_0]$ 和 $FWT[A_1]$，考虑如何合并求出 $FWT[A]$。

有以下公式：

$$FWT[A]=\operatorname{merge}(FWT[A_0],FWT[A_0]+FWT[A_1])$$

其中 $+$ 表示将序列的各位相加，$\operatorname{merge}$ 表示拼接两个序列。

这个公式的正确性大概比较显然吧（

逆变换（IFWT）

要求 $IFWT(FWT(A))=A$。把变换操作反过来即可。

$$IFWT[A]=\operatorname{merge}(IFWT[A_0],IFWT[A_1]-IFWT[A_0])$$

代码实现

两种操作只有最后一个符号不同，可以写在同一个函数中。

vector<int> FWT_or(vector<int> a, int k)
{
    if(a.size() == 1) return a;
    vector<int> a0, a1, fwta0, fwta1, ans;
    for(int i = 0; i <= a.size() / 2 - 1; i++) a0.push_back(a[i]);
    for(int i = a.size() / 2; i <= a.size() - 1; i++) a1.push_back(a[i]);
    fwta0 = FWT_or(a0, k), fwta1 = FWT_or(a1, k);
    for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back(fwta0[i]);
    for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back(fwta1[i] + k * fwta0[i]);
    return ans;
}

与

与操作与上文的或操作基本一致。

构造

我们要求：

$$C_i=\sum_{j \& k=i}A_jB_k$$

我们构造：

$$FWT[A]_i=\sum_{i|j=j}A_j$$

构造正确性的证明

$$\begin{aligned} &~~~~~FWT[A]_i \times FWT[B]_i\\ &= \sum_{i|j=j}A_j \times \sum_{i|k=k}B_k\\ &= \sum_{i|j=j}~\sum_{i|k=k} A_jB_k\\ &= \sum_{i|(k\&j)=(k\&j)} A_jB_k\\ &= \sum_{i|x=x}~\sum_{k\&j=x} A_jB_k\\ &= \sum_{i|x=x} C_x\\ &= FWT[C]_i \end{aligned}$$

变换（FWT）

$$FWT[A]=\operatorname{merge}(FWT[A_0]+FWT[A_1],FWT[A_1])$$

逆变换（IFWT）

$$IFWT[A]=\operatorname{merge}(IFWT[A_0]-IFWT[A_1],IFWT[A_1])$$

代码实现

vector<int> FWT_and(vector<int> a, int k)
{
    if(a.size() == 1) return a;
    vector<int> a0, a1, fwta0, fwta1, ans;
    for(int i = 0; i <= a.size() / 2 - 1; i++) a0.push_back(a[i]);
    for(int i = a.size() / 2; i <= a.size() - 1; i++) a1.push_back(a[i]);
    fwta0 = FWT_and(a0, k), fwta1 = FWT_and(a1, k);
    for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back(fwta0[i] + k * fwta1[i]);
    for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back(fwta1[i]);
    return ans;
}

异或

构造

我们要求：

$$C_i=\sum_{j \oplus k=i}A_jB_k$$

我们用 $x\circ y$ 表示 $\operatorname{popcount}(x\&y)\mod 2$。

我们构造

$$FWT[A]_i=\sum_{i\circ j=0}A_j - \sum_{i\circ j=1}A_j$$

构造正确性证明

可以证明：

$$(x\circ y)\oplus(x\circ z) = x \circ (y\oplus z)$$

所以：

$$\begin{aligned} &~~~~~FWT[A]_i \times FWT[B]_i\\ &= (\sum_{i\circ j=0}A_j - \sum_{i\circ j=1}A_j)(\sum_{i\circ k=0}B_k - \sum_{i\circ k=1}B_k)\\ &= (\sum_{i\circ j=0}\sum_{i\circ k=0}A_jB_k+\sum_{i\circ j=1}\sum_{i\circ k=1} A_jB_k)-(\sum_{i\circ j=0}\sum_{i\circ k=1}A_jB_k+\sum_{i\circ j=1}\sum_{i\circ k=0}A_jB_k)\\ &= \sum_{i\circ(j\oplus k)=0}A_jB_k-\sum_{i\circ(j\oplus k)=1}A_jB_k\\ &= (\sum_{i\circ x=0}\sum_{j\oplus k=x}A_jB_k)-(\sum_{i\circ x=1}\sum_{j\oplus k=x}A_jB_k)\\ &= \sum_{i\circ x=0}C_x-\sum_{i\circ x=1}C_x\\ &= FWT[C]_i \end{aligned}$$

变换（FWT）

$$FWT[A]=\operatorname{merge}(FWT[A_0]+FWT[A_1],FWT[A_0]-FWT[A_1])$$

证明：

我们再观察式子：

$$FWT[A]_i=\sum_{i\circ j=0}A_j - \sum_{i\circ j=1}A_j$$

• 对于前半段的某个位置 $i$，当前二进制位的值是 $0$。那么 $j$ 当前二进制位的取值不会影响到 $i\circ j$ 的值。所以此时 $FWT[A]_i=FWT[A_0]_i+FWT[A_1]_i$。
• 对于前半段的某个位置 $i$，当前二进制位的值是 $1$。若 $j$ 在前半段，则 $j$ 当前二进制位的值是 $0$，不会影响到 $i\circ j$ 的值。若 $j$ 在后半段，则 $j$ 当前二进制位的值是 $1$，$i\circ j$ 就会相反。所以此时 $FWT[A]_i=FWT[A_0]_i-FWT[A_1]_i$。

逆变换（IFWT）

$$IFWT[A]=\operatorname{merge}(\frac{IFWT[A_0]+IFWT[A_1]}{2},\frac{IFWT[A_0]-IFWT[A_1]}{2})$$

代码实现

vector<int> FWT_xor(vector<int> a, int k)
{
    if(a.size() == 1) return a;
    vector<int> a0, a1, fwta0, fwta1, ans;
    for(int i = 0; i <= a.size() / 2 - 1; i++) a0.push_back(a[i]);
    for(int i = a.size() / 2; i <= a.size() - 1; i++) a1.push_back(a[i]);
    fwta0 = FWT_xor(a0, k), fwta1 = FWT_xor(a1, k);
    if(k == 1)
    {
        for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back(fwta0[i] + fwta1[i]);
        for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back(fwta0[i] - fwta1[i]);
    }
    else
    {
        for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back((fwta0[i] + fwta1[i]) / 2);
        for(int i = 0; i <= a.size() / 2 - 1; i++) ans.push_back((fwta0[i] - fwta1[i]) / 2);
    }
    return ans;
}

多项式乘法逆

题目链接：洛谷 P4238

题目大意

给定 $n$ 项的多项式 $F(x)$，求出多项式 $G(x)$，使得 $F(x)\times G(x) \equiv 1 \pmod{x^n}$。即 $F(x)\times G(x)$ 只保留次数低于 $x$ 的项后的结果是 $1$。答案的系数对 $998244353$ 取模。

思路

我们假设 $n$ 是二的整数次方。我们已经递归求出了模 $x^{\frac{n}{2}}$ 的逆元 $G'$，考虑如何推出模 $x^n$ 的逆元 $G$。

$$F\times G' \equiv 1 \pmod{x^{\frac{n}{2}}}$$

$x^n$ 是 $x^{\frac{n}{2}}$ 的倍数。$F\times G$ 和 $1$ 在模 $x^n$ 下同余，则它们在模 $x^{\frac{n}{2}}$ 下也一定同余：

$$F\times G \equiv 1 \pmod{x^{\frac{n}{2}}}$$

两式相减，得

$$F\times G' - F\times G \equiv 0 \pmod{x^{\frac{n}{2}}}$$

左右同时除以 $F$，得

$$G' - G \equiv 0 \pmod{x^{\frac{n}{2}}}$$

左右平方，得

$$G'^2 - 2GG' + G^2 \equiv 0 \pmod{x^n}$$

左右同时乘 $F$，得

$$FG'^2 - 2G' + G \equiv 0 \pmod{x^n}$$

所以

$$G \equiv 2G' - FG'^2 \pmod{x^n}$$

实现

按照上面的式子直接写就行了。

因为带有取模，我们需要使用 NTT。

void solve(int n, int *a, int *b) // b 是模 x^n 下 a 的逆元
{
    if(n == 1)
    {
        b[0] = qpow(a[0], p - 2);
        return;
    }
    solve(n >> 1, a, b);
    int f[n << 1], g[n << 1], h[n << 1];
    for (int i = 0; i < (n << 1); i++) f[i] = g[i] = h[i] = 0;
    for(int i = 0; i < n; i++) f[i] = a[i], g[i] = b[i];
    NTT(n << 1, f, 1), NTT(n << 1, g, 1);
    for (int i = 0; i < (n << 1); i++) h[i] = (2 * g[i] % p - 1ll * f[i] * g[i] % p * g[i] % p + p) % p; // 这一行就是上面的式子
    NTT(n << 1, h, -1);
    int inv = qpow(n << 1, p - 2);
    for(int i = 0; i < n; i++) b[i] = 1ll * h[i] * inv % p;
}

int a[N], b[N], n;

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n;
    for(int i = 0; i < n; i++) cin >> a[i];
    int k = 1;
    while(k <= n) k *= 2;
    solve(k, a, b);
    for(int i = 0; i < n; i++) cout << b[i] << " ";
    return 0;
}

多项式除法

题目链接：洛谷 P4512