多项式

发表于 更新于 分类于 Waline：
摘要：多项式学习笔记

快速傅里叶变换（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$ 次幂下的点值，即：

暴力做是 $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(x)=A_0(x^2)+xA_1(x^2)$。

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

单位根有性质 $\omega_n^k=\omega_{2n}^{2k}$，可以推出 $\omega_n^{2k}=\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$ 的变换结果。时间复杂度 $O(n \log n)$。

1
2
3
4
5
6
7
8
9
10
11
12
13
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}$。

先写公式：

证明

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

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

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

1、$i=k$

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

2、$i \neq k$

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

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

综上，我们可以得到：

所以：

证毕。

递归实现

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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$，我们来模拟它的递归过程：

1
2
3
4
第一层：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，我们来找它们的规律。我们先写出它们的二进制：

1
000 100 010 110 001 101 011 111

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

1
000 001 010 011 100 101 110 111

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

1
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)$ 计算出。

1
2
3
4
5
6
7
8
9
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完整代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#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：

性质 3：

由费马小定理得：

所以 $\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$。

那么：

性质 4：

由等比数列求和公式得：

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

实现

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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。

快速沃尔什变换（FWT）

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

即 $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$。

或

构造

我们要求：

我们的构造方案是：

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

构造正确性的证明

变换（FWT）

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

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

有以下公式：

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

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

逆变换（IFWT）

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

代码实现

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

1
2
3
4
5
6
7
8
9
10
11
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;
}

与

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

构造

我们要求：

我们构造：

构造正确性的证明

变换（FWT）

逆变换（IFWT）

代码实现

1
2
3
4
5
6
7
8
9
10
11
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;
}

异或

构造

我们要求：

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

我们构造

构造正确性证明

可以证明：

所以：

变换（FWT）

证明：

我们再观察式子：

  • 对于前半段的某个位置 $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）

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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$。

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

两式相减，得

左右同时除以 $F$，得

左右平方，得

左右同时乘 $F$，得

所以

实现

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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