多项式的计算(上)

本文约 14300 字。

Key Words:多项式乘法,加法卷积,分治乘法,单位根,点值,插值,线性算法,快速傅里叶变换 FFT,快速数论变换 NTT,复域压缩优化,多项式泰勒展开,牛顿迭代,多项式全家桶(初等函数),分治 FFT,半在线卷积,全在线卷积,拉格朗日插值,多项式多点求值,多项式快速插值

多项式乘法

  • 多项式

    对于有限数列 ,形如 的函数称为多项式。

用 表示 的 项系数,或简写为 。

多项式乘法是多项式计算的基础,其快速算法的重要性不言而喻。接下来,我们将关注有限多项式的乘法问题。即:给出 ,,求 的各项系数。

朴素算法

根据 ,将 展开得

为了方便书写,我们认为 以及之后的项都是 ,对于 同理。

两侧同时提取 系数,得

这阐明了多项式乘法对应的数列运算。

  • 卷积

    某二元运算 的定义域和值域均为 ,数列 的“ 卷积”为一个新数列 ,满足

不难发现,多项式乘法对应系数数列的“加法卷积”。容易 暴力计算。

1
2
3
4
5
6
7
void mul(int *A, int *B, int *C, int n) {
for (int k=0; k<2*n-1; k++)
C[k] = 0;
for (int i=0; i<n; i++)
for (int j=0; j<n; j++)
C[i+j] += A[i]*B[j];
}

列竖式是手动计算加法卷积的一种简便方法,具体操作如下

分治乘法*

不妨设 为偶数,不足则末尾补 。

将 左右拆半为 ,其中

类似将 左右拆半为 ,则

(将 简写为 ,余者类似)

只需分别计算出 的各项系数,就可以用简单的加法得到答案。

一种显然思路是分别计算 ,转化为 个规模为 的子问题。复杂度递推式 ,解得 ,和暴力相同。

一种巧妙的思路是计算

则

这样可以转化为 个规模为 的子问题。复杂度递推式 ,解得 即 。

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
using Poly = vector<int>;
Poly operator + (const Poly &A, const Poly &B) {
int n=A.size();
Poly C(n);
for (int i=0; i<n; i++)
C[i] = A[i]+B[i];
return C;
}
Poly mul(Poly A, Poly B) {
int n = A.size();
if (n==1)
return (Poly){A[0]*B[0]};
if (n&1) {
A.push_back(0);
B.push_back(0);
n++; // 若 n 是奇数,则末尾补 0
}
Poly A0(n/2), A1(n/2), B0(n/2), B1(n/2);
for (int i=0; i<n/2; i++){
A0[i] = A[i];
A1[i] = A[i+n/2];
B0[i] = B[i];
B1[i] = B[i+n/2];
}
Poly Q1 = mul(A0, B0),
Q2 = mul(A1, B1),
Q3 = mul(A0+A1, B0+B1);
Poly C(n*2-1);
for (int i=0; i<n-1; i++) {
C[i] += Q1[i];
C[i+n/2] += Q3[i]-Q1[i]-Q2[i];
C[i+n] += Q2[i];
}
// C = Q1 + (Q3-Q1-Q2)*x^(n/2) + Q2*x^n
return C;
}

该示例大量使用了 ,常数开销大,仅便于理解,实用性较差。

求值-插值法

  • 点值 将多项式 视为函数,对于常数 ,记 为 在 处的点值。

  • 插值 只知道 的若干点值,求出 的系数序列。

多项式插值相当于解线性方程组(即“待定系数法”),根据线性代数知识,还原一个 次多项式( 个系数)至少需要 个不同的点值。

为了理解点值的性质,看看下面的例子:

  • 多项式乘法的检验

    给出 ,判定是否有 。

随机一个常数 ,代入算得 ,若 ,则 。

若 ,说明 。重复多次,若点值均正确,则 大概率成立。


我们已经见识了点值的妙用:通过选定具体的常数,将复杂的多项式乘法“浓缩”为数的乘法。基于这一思想,可以另辟蹊径,设计多项式乘法新的算法框架。

  • 第一步:选择足够大的 ,选定 个不同的数
  • 第二步:求出点值 与
  • 第三步:根据 得到
  • 第四步:插值求出 的系数

通过求值和插值,我们将多项式 和 转化为了 个点值,计算点值的乘法要比计算多项式的乘法要简单的多。

对于一般的 ,求值和插值较为困难。不过, 可以任我们选择,选择一组具有特殊性质的 ,可以大大简化计算。这些理想中的 ,就是单位根。

单位根

  • 复数

    令虚单位 满足 ,形如 ()的数为复数。

  • 复数的三角表示

    将复数 看做平面直角坐标系中的点 ,记 ,称 为 的模长,射线 与 轴正半轴的夹角 称作 的幅角。有 。

  • 棣莫弗定理(复数乘法的几何性质)

    两个复数相乘,模长相乘,幅角相加。

证明:设两个复数为 。

则

根据正弦、余弦的和角公式,。

(如果你是初中生,也可以用 两点间距离公式 / 相似 证明)


  • 本原单位根

    次本原单位根记为 ,它满足

    • 互不相同。
    • 。

次本原单位根可能有多个,在 中,一种经典的构造是 这是模长为 ,幅角为 ( 圆周)的复数。不难发现,乘 相当于绕着原点旋转 圆周,乘 次恰好回到原地。容易验证它满足单位根所需的性质。

  • 单位根的性质
    • 互不相同
    • (推论:)
    • (一般情况:)

根据几何性质,不难理解它们。

快速傅里叶变换 Fast Fourier Transform

我们将利用单位根补完“求值-插值法”的框架:

  • 第一步:选定 个数 (由单位根的性质,这些数互不相同)
  • 第二步:对 求出点值 与 (这一步称作 DFT)
  • 第三步:根据 得到 的点值 (这一步称作点积)
  • 第四步:插值求出 的系数 (这一步称作 IDFT,即 DFT 的逆过程)

为了便于分治,考虑 为 的整数次幂的情况(不足则末尾补 )。

DFT

给出 ,下面我们将对 求出 。

考虑分治,将 奇偶拆半为 其中 含偶数项系数, 含奇数项系数。它们可以表示出 。

假设我们已经知道对 求解子问题的答案 为了求原问题的答案,将 代入。 此处用到的点值我们都已经取得,故可以用子问题的答案求出原问题的答案。

使用分治,复杂度递推式为 ,解得 。

为什么是奇偶拆半而不是左右拆半?这本质上是整除 转化为子问题,然后讨论余数合并。不一定要分成 个子问题,可以分成 个或更多。

DFT 的线性性与 IDFT

接下来先引入线性算法的概念,再说明 DFT 是线性算法,最后对线性算法求逆得到 IDFT。

  • 线性算法

    线性算法的输入为数列 ,输出为数列 。算法有一个内秉的系数矩阵 ,算法负责计算

记 是一个线性算法, 是一个数列,以 为输入计算 的输出记作 。

根据线性代数知识,线性算法相当于“向量乘矩阵得到结果向量”,记作 。如果要得出逆算法,只需要求出 的逆矩阵 ,然后就有 了。

  • 线性算法的线性性

    对于线性算法

    • (其中 为常数)

在具体应用中,常利用 DFT(IDFT)的线性性减小常数。

观察 不难发现 DFT 是线性算法,其转移矩阵 ,即

现在让我们暂时忘掉多项式,将 DFT 视作“输入一个数列,输出一个数列的算法黑箱”。我们希望求出其逆算法,而这归根结底是寻找逆矩阵,有结论

证明:

记 ,,只需证明两矩阵相乘为单位阵 。

记 ,有 这是等比数列求和。

若 ,每一项都是 ,和是 ,。

若 ,运用等比数列求和公式,和是 ,其中 ,故 。

综上 ,。


我们已经得到了可行的 IDFT 算法,如何快速计算?注意到 和 非常相似,IDFT 可以转化为 DFT。

欲对数列 计算 ,有

构造多项式 ,则 。只需求解 ,这和 DFT 非常相似。

  • 实现方法一:将单位根的定义修改为 ,也符合单位根所需的所有性质。在新定义下求解 DFT 再乘 ,就相当于完成了 IDFT。

  • 实现方法二:注意到 ,先做一次 DFT,将 位置的答案翻转,再乘 就可以达成 IDFT。

至此,我们成功补完了“求值-插值法”的大框架,实现了 的多项式乘法。

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
const double Pi = acos(-1);
using Complex = std::complex<double>;
using Poly = std::vector<Complex>;
// 使用 c++ 自带的 complex 类
Poly DFT(const Poly &F) {
if (F.size()==1)
return F;
int n = F.size();
Poly F0(n/2), F1(n/2);
for (int i=0; i<n/2; i++){
F0[i] = F[i*2];
F1[i] = F[i*2+1];
}
Poly S(n), S0 = DFT(F0), S1 = DFT(F1);
Complex w(cos(2*Pi/n), sin(2*Pi/n)),
buf(1,0);
for (int k=0; k<n; k++){
S[k] = S0[k%(n/2)] + buf*S1[k%(n/2)];
buf = buf*w; // 此处 buf = w_n^k
}
return S;
}
Poly IDFT(const Poly &F) {
Poly S = DFT(F);
int n = S.size();
reverse(&S[1], &S[n]);
// 采用方法二,将 S[1] ~ S[n-1] 翻转
for (int i=0; i<n; i++)
S[i] /= n;
return S;
}
Poly mul(Poly A, Poly B) {
int M = A.size()+B.size()-1, N;
for (N=1; N<M; N<<=1);
A.resize(N);
B.resize(N);
// 将 N 补至 2 的整幂
A = DFT(A);
B = DFT(B);
for (int i=0; i<N; i++)
A[i] = A[i]*B[i];
A = IDFT(A);
return A;
}

DFT 的常数优化

  • 单位根复用

注意到 ,讨论余数将取模展开,对于 在计算两者时, 可以共用,这样能节省半数的乘法。

此外,我们用到的单位根只有 共 个。但我们中途用累乘法计算这些单位根,实际上使用了 次乘法。预处理这些单位根,不仅能减少时间开销,还能降低浮点数精度损失。

  • 蝴蝶变换和按层迭代

观察分治树,整个计算过程可以分为两步:

  • 第一步:反复执行“奇偶拆半”,将 的系数分到叶节点上。
  • 第二步:从下往上计算,利用两个儿子的信息计算出父亲的信息。

image-20250324151237018

第一步的流程如上图。记 , 为将 的二进制位反转的结果( 视为一个 位二进制数,含前导 )。可以发现,系数 最终到达了位置 。掌握了这个规律,我们可以用 的简单交换(而非 的反复奇偶拆半)来完成第一步。

对于第二步,原有的递归算法涉及到大量 的传递,开销较大。可以将递归分治改为逐层迭代执行。具体地,将一层所有分治树节点的答案从左到右存储在一个数组中,每次向上计算一层。最开始叶子层的结果由蝴蝶变换给出。

优化后的代码如下。该实现代码简短,便于记忆,为大部分选手所采用。

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
const int MaxN = ______;
using Complex = std::complex<double>;
Complex w[MaxN];
int rev[MaxN];
void Init(int N) {
for (int n=2; n<=N; n*=2)
for (int k=0; k<n/2; k++)
w[n/2+k]=(Complex){cos(2*Pi*k/n),sin(2*Pi*k/n)};
// 预处理单位根,将 w_n^k 存储在 w[n/2+k],巧妙地避免了冲突
for (int i=1; i<N; i++)
rev[i] = (rev[i>>1]>>1) | ((i&1) ? N/2 : 0);
// 递推求出 rev 数组
}
void DFT(Complex *F, int N) {
for (int i=0; i<N; i++)
if (i<rev[i])
swap(F[i], F[rev[i]]);
// 应用蝴蝶变换,预处理叶子层的元素
for(int n=2; n<=N; n<<=1) {
// n 为当前层的子问题规模
for(int l=0; l<N; l+=n) {
// l 为当前子问题的左侧起始位置
// F[l ~ l+n/2-1], F[l+n/2 ~ l+n-1] 是左右儿子的答案
for(int k=0; k<n/2; k++) {
Complex sav = w[n/2+k]*F[l+k+n/2];
F[l+k+n/2] = F[l+k]-sav;
F[l+k] = F[l+k]+sav; // 注意计算顺序
}
}
}
}

复域压缩优化

利用 DFT 算法的性质,可以节约计算 DFT 的次数。

  • 三次变两次

    给出 ,欲计算 。

    构造复多项式 ,则 ,只需计算 ,取出其虚部就能得到 的系数。而计算 只需要一次 DFT 和一次 IDFT。(朴素算法需要两次 DFT、一次 IDFT)

  • 两次变一次*

    • 复数的共轭

      称复数 与 (实部相同,虚部互为相反数)互为共轭(镜像对称)。记 为复数 的共轭。

    • 共轭的性质

      复数的运算具有镜像对称性,对于复数

    这种对称性允许我们将整个线性算法共轭。

    • 共轭算法的性质

      对任意线性算法 ,记 为它的共轭算法,两者的转移矩阵互为共轭。则

    将线性算法展开易证。


    考虑如何计算 ,易知 ,注意到 ,根据前文对 IDFT 的讨论,对 做 DFT 再将 位翻转相就当于 。

    现给出 ,欲计算 。

    构造一对共轭多项式 先用计算 ,然后 因此,将 的 位翻转,再取共轭就得到 。

    根据线性性, 相加减(并提取实部、虚部)就能分别得到 ,仅需一次 DFT。

    若要同时计算两个 IDFT,方法类似。

快速数论变换

许多计数题都要求答案对一个较大的整数 取模,如果在模 意义下计算多项式乘法,可以先视作值域为 的普通多项式乘法,再对结果模 。

答案系数可达到 级别,若取 ,则达到 ,超出 的精度范围,直接用 FFT 计算会产生精度问题。

如果能在同余系中找到一个数,满足单位根的性质,就能绕开复数,直接在模意义下计算多项式乘法了。

在一些非常特殊的模数下,这的确可以实现。例如质数 ,其原根 ,若 是 的幂,构造 ,其阶为 ,完全符合单位根所需的性质。

是多项式计数题中最为常见的模数。

任意模数多项式乘法:拆系数 FFT

不利用模数的特殊性质。考虑直接先相乘后取模,这需要设计一个精度足够高的多项式乘法。

题目中的系数 一般在 以内,将其拆成 的形式,其中 。

将多项式拆成 ,其中 的系数都不超过 ,若其中两个相乘,结果的值域仅仅是 级别, 可以承受。

我们要计算的是

暴力计算需要 次 DFT(用于 )和 次 IDFT(用于 ),共 次。可以用“两次变一次”优化到 次。此外,若不想借助“两次变一次”,也可以仿照分治乘法优化成 次乘法,再用“三次变两次”,得到共 次(或 次) DFT/IDFT。

还有一种解决方案是三模数 NTT,需要中国剩余定理合并答案,较为繁琐,故不介绍。


  • 例题 6.1.1. 差卷积

    给出两个数组 ,计算两者的差卷积 ,其中

Solution:将 翻转得到 ,使得 。则

第三个等号改为枚举 。

转化为加法卷积,将 卷积后取出 项即可。复杂度 。



形式幂级数

本章名为“多项式计数”,是一种通俗的叫法,实际上我们分析的主要是形式幂级数,在实际计算时归为多项式。

形式幂级数的引入

在这一节中,我们将有限项多项式扩展到无限项的情形。首先定义无穷数列的运算。

  • 无穷数列的运算

    对于无穷数列 ,它们的加法、数乘和加法卷积定义为:

以这些运算为基础,我们可以将(有限项)多项式扩展为(无穷项的)形式幂级数。

  • 形式幂级数

    对于无穷数列 ,形如 的式子称为形式幂级数。

  • 形式幂级数的运算

    对于形式幂级数 ,,它们的加法、数乘和乘法定义为:

这三个运算的定义和多项式一致。当无穷序列只有有限前缀非零时,形式幂级数退化为多项式。为了分析形式幂级数,可以从其有限前缀入手,这需要引入截断的概念。

  • 截断

    表示保留 的 次项得到的多项式。

两个形式幂级数 相等,当且仅当它们每一项都相等,即对于任意大的 ,。

  • 截断的性质

    对于形式幂级数 ,有

这两项性质容易验证。它们说明,如果我们只对形式幂级数的前 项感兴趣,计算时只需保留前 项。

形式幂级数和幂级数的区别*

  • 幂级数

    对于无穷数列 ,形如 的函数称为幂级数。

例如我们熟悉的几何级数 可以证明该无穷求和在 时收敛,且

“形式幂级数”和“幂级数”的分歧在于,后者关注无穷求和的取值,而前者关注系数数列。例如对于 若将它们视作幂级数,都只在 时才收敛,即 未定义 这毫无意义。然而,若将它们视作形式幂级数,那么它们是不同的,它们分别代表数列 换句话说,幂级数中的 代表“某个未知的实数”,而形式幂级数的 只是一个记号, 用于标记数列的第 项,此外别无他用。

未来的章节将关注形式幂级数,因为我们关注系数数列而非无穷级数的和。我们需要关注系数的收敛性,而非整个无穷求和的的收敛性。

例如对于 若将其视作形式幂级数,考察常数项 发现其不收敛。但是,若将其视为无穷求和,则有 这再次印证了两者的差异。

  • 表:形式幂级数与幂级数的对比
形式幂级数 幂级数
本质 无穷数列 函数(无穷求和)
的含义 用于标记数列项 某个未知实数
收敛性要求 系数收敛 无穷求和收敛

形式幂级数的初等函数

从形式幂级数的乘法出发,可以定义一系列初等函数。

  • 形式幂级数的幂
    • 乘法逆:若 ,称 互为乘法逆元。
    • 有理次幂:若 ,称 。

先前我们总是将形式幂级数写成逐项求和的形式,例如 可以发现 因此可以记 其中 被称作“封闭形式”。

“封闭形式”指不含不定求和或不定乘积等的简单形式。它是一种约定俗成的概念,例如“”暗含了不定乘积,仍被视为一种可接受的封闭形式。


  • 形式幂级数的复合

    对于形式幂级数 ,记它们的复合为

两个形式幂级数的复合会产生无穷求和,这可能导致收敛性问题。

  • 复合收敛的充分条件

    对于形式幂级数 ,若 ,则 一定收敛。

证明:由于 , 的最低项至少是 。 后者是有限求和,必然收敛。


  • 形式幂级数的指对函数

    定义形式幂级数 然后用复合来定义指数和对数。对于形式幂级数 ,,定义

和 的定义和实函数 , 的泰勒展开式相同。它们也拥有指数和对数的运算性质,例如


  • 形式幂级数的导数和积分

    对于形式幂级数 ,定义

    • 形式导数

    • 形式不定积分 其中 为任意常数。

可以验证上述记号符合我们想要的微积分性质,例如


接下来,我们将给出计算各类初等函数的算法。

乘法逆

给出 , 满足 ,求 。

  • 递推法

对比系数知 。考虑递推,设 且 已经求出,欲求 。

边界是 ,可以 递推。

  • 倍增法

假设已经计算出 ,考虑如何计算 。 于是,可以由 计算出 ,每次使次数加倍。复杂度递推式 解得复杂度为 。后文将沿用这一复杂度分析。

对数函数

给出 的前 项,满足 ,求 。

设 ,两侧求导得 积分得 不定积分不确定常数项,此时根据定义易知 。

牛顿迭代法

传统意义上的牛顿迭代法可以求 上函数方程的数值解。

对于函数 ,欲求方程 的数值解。先取猜测值 ,构造如下递推: 该构造有一个简单的几何解释:在 处作切线,取切线与坐标轴的交点作为新的估计值。下图展示了两步牛顿迭代的效果。

随着迭代的进行, 会逐渐接近方程的解。可以证明牛顿迭代在某些前提下拥有平方收敛的性能,粗略地说,每迭代一次,结果的有效数位增加一倍。


记形式幂级数的集合为 , 上的牛顿迭代可以扩展到 上。

  • 形式幂级数上的函数

    可以令函数 为定义在形式幂级数上的函数,其定义域和陪域都为 。

例如 ,根据定义,对于多项式(或形式幂级数) ,。

  • 形式幂级数方程

    给出函数 ,我们的目标是解方程 。它能确定形式幂级数 ,在实践中,我们感兴趣的是它的截断 。

例如 ,那么 的解 就是 的乘法逆元。

一般地,如果 中只含有熟悉的初等函数,可以用形式幂级数的复合来定义它。例如 可以视为二元形式幂级数,有

  • 形式幂级数(复合)方程的截断

    可以用复合定义。对于方程 ,记 ,有

利用截断的性质易证。我们称 为方程 在 精度下的解。


  • 上的泰勒展开

    若 可以用复合定义,定义 的导数为二元形式幂级数的偏导 对于形式幂级数 ,若 收敛,则

证明:先证 上的麦克劳林公式,即 的情形

记 计算得 因此 麦克劳林公式得证。令 (可证其收敛),对其应用麦克劳林公式: 再令 得


理论准备完毕,下面介绍如何求方程 在 精度下的解。

  • 上的牛顿迭代

    欲求 在 下的解,假设已知 ,则

证明:将 在 处泰勒展开,再代入 得

由 , 的最低次项至少是 ,可知 等的最低次项至少是 。对上式模 ,这些项都变成 ,故 整理即得上文的倍增公式。


利用该公式,每次可以将解的次数倍增。其形式和 上的牛顿迭代公式一致,可以借此记忆。 上的“平方收敛”对应多项式的“次数倍增”。

在计算时需要注意,由于 的前 项是 ,分母 的精度只需达到 就足够了。

利用这个结果,再次探究多项式求逆的倍增算法。构造 ,欲解 ,假设已知 ,套用牛顿迭代公式 注意 的精度只需达到 ,用 代替即可,最终得到的倍增公式与前文相同。

多项式开平方

给出 , 满足 ,求 。

构造 ,则 ,。

假设已求出 满足 ,套用倍增公式整理得

复杂度 。

多项式指数函数

给出多项式 ,记 ,求 。

递推法:求导得 ,提取 系数得 可以 求出末项,复杂度 。

将 分离可以得到多项式 的递推方法

倍增法:构造 ,我们要求 的解。

假设已求出 满足 ,套用牛顿迭代公式整理得 复杂度 。

多项式 k 次幂

给出 和 ,求 。

使用快速幂计算,复杂度 。

若 ,则 ,复杂度 。

若 ,找出 的最低次项 ,计算 的幂(此时满足 ),再乘以 即可。


当 项数较少时,可以用递推法快速计算。

记 ,两侧求导可得 ,两侧同乘 得 ,提取 系数得

其中 只需要枚举 中有值的项,设有 个,复杂度 。

注:实际上这套方法是“整式递推”,更多知识可见后文《递推关系》一节。

多项式带余除法

  • 多项式取模

    给出 ,其中 是首一多项式。若多项式 满足 ,且 的次数小于 ,则称 为 整除 的余数,记为 。

  • 多项式取模的性质

    对于多项式 ( 是首一多项式),有

接下来考虑如何计算多项式取模:给出多项式 ,求 。

如上设 ,考虑用某种方式把余多项式 的干扰去除。记 次多项式 的系数翻转为 。记 的次数为 ,将 的次数视作 ,将 系数翻转得 对 取模得 这样就去掉了余数 。直接多项式求逆即可得到 ,进而由 得到 。复杂度 。

当 项数较少时,可以用递推法(模拟初中学过的竖式大除法)快速计算。


分治 FFT

多项式的连乘积

给出 个多项式 ,其次数总和为 ,求 。

采用分治法,记 ,则有转移

分治树有 层,总复杂度 。

分式求和

给出 个分式 , 的次数总和为 ,求

的前 项系数。

分别维护分子分母,,按照这个规则分治合并。

算得 ,其中 的次数为 。接下来多项式求逆即可。

半在线卷积

  • ,初始时给出 ,计算出 后,才能获得 的值。求 的前 项,该问题称作“半在线卷积”。

类似 DP 优化,使用“中序遍历的 CDQ 分治”,对于当前区间 ,记 。考虑跨越隔板的贡献,即 对 的贡献。对于 有

此时已递归左区间, 已计算完毕,故已得到 ,上式可卷积计算。

到达叶节点 时,则计算 。总复杂度 。


可以分多叉来降低复杂度,假设将当前区间均分成 段,会产生 次跨越隔板的贡献,但利用点值的线性性,只需要计算 次 DFT/IDFT,配合 次点乘。总复杂度 ,取 得复杂度为 。

一个经典的应用是多项式 。记 ,将前文递推公式变形得 这类似于半在线卷积。于是可以 或 求解多项式 。相比于牛顿迭代法,此法的思路和代码实现都更简洁,在效率上也不差。

全在线卷积

  • ,在你正确地计算出 后,才能获得 的值。求 的前 项系数。

为了便于理解,想象一个阶梯状的矩阵,位置 代表乘积 , 是某一斜行的和。

的具体流程如下:

  • 记 ,调用 。
  • 处理 的 对 的贡献,分两类讨论:
    • 若 ,计算
    • 若 ,计算
  • 调用 。

例如对于 ,三个步骤计算的贡献如图所示,分别对应灰色、蓝色和红色:

当到达叶子 :

  • 若 ,。
  • 若 ,。
  • 已得到正确的 ,可取得新的 。

可以发现,算法不会用到未知的 ,符合全在线卷积的要求。

调用 就能完成计算。复杂度 得 。


还有一种等价算法。一个规模为 的全在线卷积可以按顺序拆分成:

  • 一个规模为 的全在线卷积。(灰)
  • 一个规模为 的普通卷积。(蓝)
  • 两个并行的,规模为 的半在线卷积。(红)

复杂度递推式 得 。

对于上述两种算法,多叉平衡均可以让复杂度更优。


求值与插值

拉格朗日插值

给出 次多项式 的 个点值 ( 互不相同),求 的各项系数。

这相当于解线性方程组

首先证明解唯一。根据线性代数知识,相当于证明系数矩阵 满秩,它是一类特殊矩阵:

  • 范德蒙德矩阵

    称各行为几何级数的矩阵 为 对应的范德蒙德矩阵。

证 满秩即证 。

  • 范德蒙德行列式

    对于 对应的范德蒙德矩阵 ,有

可以用归纳法证明。

由 互不相同,易知 ,因此解是唯一的。


接下来尝试构造出解。先解决 的情形,不难发现

对于一般情况,首先对 构造 次多项式 使得 ,,则 先构造 ,则 ,且 ,则 综上有插值公式 将 展开得

  • 拉格朗日插值公式

    对于互不相同的



多项式多点求值

给出多项式 与常数 ,求各个 的值。

  • 引理

证明: 是一个常数,记为 ,则 。代入 得 。


记 ,我们的目标是对每个 求出 。

根据多项式取模的性质,记 ,对 计算 时,可以先对 模 ,这不会影响答案。

于是有如下分治算法:

  • 对于分治区间 ,记 。假设目前已经有 。
  • 计算 ,分治处理区间 。
  • 计算 ,分治处理区间 。
  • 到达叶子 时,能获得 。

其中 可以用一次自下而上的分治 FFT 预处理。在子问题分治区间 时, 的次数被 削减到了 ,故总复杂度 。

对 进行多点求值的流程如下图:

下文《转置原理》一节会介绍一个常数更小、更易实现的多点求值算法。

多项式快速插值

写出拉格朗日插值公式

先对每个 求出 。

记 ,可以分治 FFT 求出。对 代入 就得到 。直接代入会使分子分母均为 ,需要使用 形洛必达法则 或利用导数的乘法法则,也可得到 于是,对 多点求值可得到各个 。

记 ,问题转化为 这可以分治 FFT 求解。对当前区间 ,记 考虑“缺口”在左侧还是右侧,有转移 复杂度 。