Halberd Cease's Blog
多项式学习笔记
Halberd Cease
  2023 年 11 月 24 日 14:58   2024 年 07 月 17 日 21:45    

插播：鰰的多项式学习笔记：多项式全家桶，FFT (快速傅里叶变换)

多项式乘法

两个多项式的乘积被定义为：

其中 是 和 的卷积。

朴素计算的时间复杂度是 的。

多项式的点值表示法

给定一个不超过 次的多项式 以及 个不同的点 ，令 ，则这 组 唯一确定了这个多项式 。

当使用点值表示法的时候，多项式相乘可以直接相乘而不是用朴素形式卷积，时间复杂度 。

问题：如何快速在系数表达和点值表达之间切换？

题外话：拉格朗日插值

有：

其中 为点值表示法下的点的坐标，这个公式可以在 的时间内将多项式的点值表示法转换成系数表示法。它的更加广泛的应用是给定平面上 个点，求一条过这 个点的 次直线（对其的横坐标没有特定的要求）。

​ 此处的 次，代表最高次为 次，不保证不会退化成更低次的多项式。

但这还是太慢了，和直接卷积的复杂度一样，效率不高。

考虑：特殊的点值。

单位根

Definition： 满足 的 被称作 次单位根。

显然， 是这个东西的一个解，于是 是 次单位根，如果 是偶数的话， 也是 次单位根。

想要寻找其它单位根，我们需要引入复数。

复数

Definition： 形如 的数，其中 是虚数单位，满足 ，就是复数。

这个东西长得很像向量，于是你可以把这个东西抽象到二维平面上（于是这个平面就被叫做复平面）。横轴表示实数，也就是一个数轴，在原点垂直这个数轴再作一条数轴，为纵轴，表示（纯）虚数。横轴的单位是 ，纵轴的单位是 ，于是在这个平面坐标系上的点就构成了整个复数域。

复数的运算

加法： 定义，那么 ；

减法： 同理；

乘法： 模长相乘，俯角相加（具体不想写了）。

单位根

再回来看单位根，有式子 ，那么 ，在复数域上，这个东西就表示 个 相乘，结果的模长为 ，我们又知道 其实是相等的而且大于 （显然），所以 的模长就应该是 ，也就是说，在复平面上，这个数到原点的距离为 。

所以，单位根会落在以原点为圆心，半径为 的圆上。

再考虑乘法中“俯角相加”这一条，这里的俯角指的是复平面中表示这个复数的向量，由实轴正半轴逆时针旋转得到的角。

还是原来那个式子，由于每一个 的弧度是相等的，而相加为 ，故 的弧度应该是 。

​ 其实相加不一定是 ， 乘起来可能转了很多圈然后回到 ，所以 的弧度实际上应该是 ，如果 的话就是和前面重复的，故不讨论。

本原单位根

Definition： 设 是 次单位根，如果 恰好生成了所有的 次单位根，则称 是本原单位根。

在复数域 上， 是一个本原单位根。

离散傅里叶变换（DFT）

设 是长度为 的数列，对 ，令

则称 为 的离散傅里叶变换，暂时记作 。

离散傅里叶变换的逆变换（IDFT）

对于长度为 的序列 ，有 如下：

其中 。

重要的，两者互逆，所以已知多项式在单位根处的点值， 可以求出多项式的各项系数，这个过程也可以看作插值。

但是，这个东西的复杂度也是 的，所以我们还需要优化。

循环卷积

对于两个长度为 的序列，定义其循环卷积为

记作 。

注意：循环卷积的长度相较原序列不变，且两个原始序列的长度需要相等。

关于循环卷积和 ，有如下的定理：

其中 表示逐点乘积。

快速傅里叶变换

是一种快速计算 的方法，时间复杂度为 。

当 为 的幂次的时候，使用 算法实现 。

实际上，前面说过，最高次项系数为 ，即多项式退化的情况不影响我们的卷积过程，所以在 OI 中任何的多项式都可以表示成项数 ， 的形式，因此都可以使用 算法去求 。

补充：单位根的一些性质

考虑 和 ，有：

蝴蝶操作

设 。考虑将 按次数的奇偶性分类：

定义

那么 和 都是次数不超过 的多项式，并且有

结合单位根的性质，对于 ，有

通过蝴蝶操作，如果我们得到了 在点 处的点值，就可以在 的时间内计算出 在 处的点值。

计算 的点值的过程可以递归分治进行。

而这个算法的时间复杂度是 的。

至此，我们已经获得了可以快速将系数表达方法的多项式转换为点值表达方法多项式的方法。

但是，最终我们还是要把计算过后的点值表达式转换回系数表达式，如何快速进行 ？

对比 和 的表达式，可以发现，只需要将 过程中的 全部替换成 ，最后再乘上 即可。

令 表示将 二进制顺序反转后得到的数字，令

那么每次需要对 进行的蝴蝶操作在 中变成了对两个相邻序列的操作。

得到 后，先对其中长度为 的相邻的点值蝴蝶操作，然后倍增地去想，对长度为 的相邻的点值蝴蝶操作，因为 是 的次幂，所以恰好能在 层操作后全部操作完，最后一次也就是最终的合操作。

Code

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
61
62
63
64
65
66
67
68
69
70
71
72
73
#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
int n, r[8000010];
struct C // 复数类
{
    double r, i;
    C() { r = i = 0; }
    C(double x, double y)
    {
        r = x;
        i = y;
    }
    C operator+(C &x) { return C(r + x.r, i + x.i); }
    C operator-(C &x) { return C(r - x.r, i - x.i); }
    C operator*(C &x) { return C(r * x.r - i * x.i, r * x.i + i * x.r); }
    void operator+=(C &x)
    {
        r += x.r;
        i += x.i;
    }
    void operator*=(C &x)
    {
        double t = r;
        r = r * x.r - i * x.i;
        i = t * x.i + i * x.r;
    }
} f[8000010], g[8000010];
void FFT(C *a, int op) // FFT
{
    C W, w, t, *a0, *a1;
    int i, j, k;
    for (i = 0; i < n; i++) // 位逆序置换
        if (i < r[i])
        {
            t = a[i];
            a[i] = a[r[i]];
            a[r[i]] = t;
        }
    for (i = 1; i < n; i <<= 1) // 循环 log n 层进行蝴蝶操作
        for (W = C(cos(PI / i), sin(PI / i) * op), j = 0; j < n; j += i << 1)
            // 这里的 *op 是运算和逆运算，上文提到过 DFT 和 IDFT 的
            // 差别其实在于乘的是 \omega_n 还是 \omega_n^{-1}，
            // 最后的乘上 1/n 在输出的时候解决，于是就可以一个函数实现
            // DFT和IDFT。
            for (w = C(1, 0), a1 = i + (a0 = a + j), k = 0; k < i; k++, a0++, a1++, w *= W)
            {
                t = *a1 * w;
                *a1 = *a0 - t;
                *a0 += t;
            }
}
int main()
{
    ios::sync_with_stdio(0);
    int m, i, l = 0;
    cin >> n >> m;
    for (i = 0; i <= n; i++)
        cin >> f[i].r;
    for (i = 0; i <= m; i++)
        cin >> g[i].r;
    for (m += n, n = 1; n <= m; n <<= 1, l++)
        ;
    for (i = 0; i < n; i++)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); // 位逆序置换的预处理
    FFT(f, 1);
    FFT(g, 1);
    for (i = 0; i < n; i++)
        f[i] *= g[i]; // 点值相乘
    FFT(f, -1);
    for (i = 0; i <= m; i++)
        printf("%.0f ", fabs(f[i].r) / n);
}

快速数论变换

和 都是在 中进行的过程，很多时候，多项式计算都是对整数操作，并且经常会对某一模数 取模。

注意单位根在 中起到了重要的作用，我们考虑在模素数意义下是否存在和单位根性质类似的元素。

设 是一素数，由费马小定理，对于任意 有

原根

Definition： 称为模 的原根，当且仅当 在模 意义下互不相同。

可以证明，模 意义下的原根总是存在的。

原根的性质和本原单位根非常类似，换句话说，在 意义下， 可以被看做一个 次本原单位根。

设 满足 。令 ，那么有

并且 互不相同。

于是 在模 意义下便有 次本原单位根的性质。用它类似地定义 ，这被称为快速数论变换（）。

的优势：不丢失精度；

的劣势：不能计算实数卷积，模数受限。

Code

不想写了，于是借鉴货车的，原文在此。

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
#include<bits/stdc++.h>
#define int long long
using namespace std;
int a[4000010],b[4000010];
int n,m;
const int mod=998244353,g=3,inv=332748118;
int qpow(int a,int b)
{
    int res=1;
    while(b)
    {
        if(b&1)res=res*a%mod;
        a=a*a%mod;
        b/=2;
    }
    return res;
}
void NTT(int *u,int len,int sign)
{
    if(len==1)return;
    int mid=len/2;
    int a1[mid+1],a2[mid+1];
    for(int i=0;i<=len;i+=2)
    {
        a1[i/2]=u[i];
        a2[i/2]=u[i+1];
    }
    NTT(a1,mid,sign);
    NTT(a2,mid,sign);
    int w=1;
    int v;
    if(sign==1)v=qpow(g,(mod-1)/len);
    else v=qpow(inv,(mod-1)/len);
    for(int i=0;i<mid;i++,w=w*v%mod)
    {
        u[i]=(a1[i]+w*a2[i]%mod)%mod;
        u[i+mid]=(a1[i]-w*a2[i]%mod+mod)%mod;
    }
}
signed main()
{
    cin>>n>>m;
    for(int i=0;i<=n;i++)
    cin>>a[i];
    for(int i=0;i<=m;i++)
    cin>>b[i];
    int base=ceil(log2(n+m+1));
    int lenth=1<<base;
    NTT(a,lenth,1);
    NTT(b,lenth,1);
    for(int i=0;i<=lenth;i++)
    a[i]=a[i]*b[i]%mod;
    NTT(a,lenth,-1);
    for(int i=0;i<=n+m;i++)
    cout<<(a[i]*qpow(lenth,mod-2))%mod<<' ';
    return 0;
}

形式幂级数

不会。

环 上的形式幂级数如

其中 。 称为这个形式幂级数的自由元。

多项式是仅有有限项非零的形式幂级数，因此，形式幂级数可以看作对多项式的推广。

一般形式幂级数的 不带入具体数值。

可以看出，形式幂级数 对应一个无穷序列 。

称 为数列 的生成函数。

在实际运算中，我们时常只关心形式幂级数的前有限项，并且希望形式幂级数在这种意义下参与运算。如果只关心前 项，我们用 表示。

设 是形式幂级数。如果 ，那么存在形式幂级数 使得 ，这时，称 是 的逆元，记作 。

生成函数

生成函数是一种用来描述特定组合对象的形式多项式，组合对象满足特定性质，如
大小，长度为 等。每一类组合对象都被定义了特定性质 ，其中满
足 的组合对象数是有限的。
组合对象通常可以分为无标号和有标号两类，其中无标号的组合对象可以直接拼接
而不考虑内部顺序，如序列等。
而有标号的组合对象的拼接则需要考虑内部顺序，如森林等。

普通生成函数（OGF）

也叫无标号生成函数，用来描述数列 的普通生成函数定义为形式幂级数

其中 是一类无标号对象， 表示满足 的组合对象数。

指数型生成函数（EGF）

有标号生成函数，数列 的指数生成函数定义为形式幂级数

形式幂级数在处理有标号组合对象时更为便捷。

考虑有标号的组合对象的拼接，将对象 拼接为 ，其中 ，由于其相对的顺序有变化，就会有

种方法，因此，若将两类组合对象 拼接为 ，则有

那么 。

完结。

附录

  1. 欧拉公式：;

    特殊形式：.

  2. .

多项式学习笔记
2023/11/24/poly/
作者
Halberd Cease
发布于
2023-11-24 14:58
许可
  1. 多项式乘法
  2. 多项式的点值表示法
    1. 题外话：拉格朗日插值
  3. 单位根
  4. 复数
    1. 复数的运算
  5. 单位根
    1. 本原单位根
  6. 离散傅里叶变换（DFT）
    1. 离散傅里叶变换的逆变换（IDFT）
  7. 循环卷积
  8. 快速傅里叶变换
    1. 补充：单位根的一些性质
    2. 蝴蝶操作
    3. Code
  9. 快速数论变换
    1. 原根
    2. Code
  10. 形式幂级数
  11. 生成函数
    1. 普通生成函数（OGF）
    2. 指数型生成函数（EGF）
  12. 附录
  1. 多项式乘法
  2. 多项式的点值表示法
    1. 题外话：拉格朗日插值
  3. 单位根
  4. 复数
    1. 复数的运算
  5. 单位根
    1. 本原单位根
  6. 离散傅里叶变换（DFT）
    1. 离散傅里叶变换的逆变换（IDFT）
  7. 循环卷积
  8. 快速傅里叶变换
    1. 补充：单位根的一些性质
    2. 蝴蝶操作
    3. Code
  9. 快速数论变换
    1. 原根
    2. Code
  10. 形式幂级数
  11. 生成函数
    1. 普通生成函数（OGF）
    2. 指数型生成函数（EGF）
  12. 附录