多项式 I:拉格朗日插值与快速傅里叶变换

1. 复数和单位根

前置知识:弧度制,三角函数。

1.1 复数的引入

跳出实数域 (mathbb R),我们定义 (i ^ 2 = -1),即 (i = sqrt {-1}),并在此基础上定义 复数 (a + bi),其中将 (bneq 0) 的称为 虚数。复数域记为 (mathbb C)

像这种从 (a) 变成 (a + bx)扩域 操作并不少见,例如初中学习 “平方根” 时,经常用 (a + bsqrt x(x > 0)) 表示一个数。这类数的加减乘都是容易的,除法即考虑平方差公式 ((c + dsqrt x)(c - dsqrt x) = c ^ 2 - d ^ 2x),因此 (frac {a + bsqrt x} {c + dsqrt x} = frac {(ac - bdx) + (bc - ad)sqrt x} {c ^ 2 - d ^ 2 x})

(x) 替换为 (-1),复数四则运算可由实数四则运算结合 (i) 的定义直接推广得到。

  • 加法:((a + bi) + (c + di) = (a + c) + (b + d)i)
  • 减法:((a + bi) - (c + di) = (a - c) + (b - d)i)
  • 乘法:((a + bi)(c + di) = (ac - bd) + (ad + bc)i)
  • 除法:(frac {a + bi} {c + di} = frac {(a + bi)(c - di)} {(c + di)(c - di)} = frac {ac + bd} {c ^ 2 + d ^ 2} + frac {bc - ad} {c ^ 2 + d ^ 2} i)

1.2 复平面与棣莫弗定理

描述一个复数 (a + bi) 需要两个值 (a)(b),其中 (a) 表示 实部(b) 表示 虚部。这启发我们将其放在平面直角坐标系上描述,称为 复平面。其在复平面上的坐标为 ((a, b)),实部 (a) 为横坐标,虚部 (b) 为纵坐标。

一个复数唯一对应一个平面向量,因为它们都可以用有序实数对描述。将向量起点平移至原点,则其终点指向与其对应的复数。考虑平面向量的一些特征,如其长度与所在直线斜率。将这些概念应用在复数上,我们得到如下定义:

  • 定义复数 (z = a + bi)(|z| = sqrt {a ^ 2 + b ^ 2})
  • 定义复数 (z = a + bi)辐角(operatorname{Arg} z = theta),其中 (tan theta = frac b a)。满足 (-pi 的 (theta) 称为 辐角主值,记作 (operatorname{arg} z),即 (operatorname{arg} z = arctan frac b a)
  • 辐角确定了 (z) 所在直线,模确定了 (z) 在直线上的长度。对比实部和虚部,模和辐角主值以另一种有序实数对的形式描述复数。

根据 (z = a + bi) 的模 (r) 和辐角 (theta),可知 (z) 的实部 (a = rcos theta),虚部 (b = r sin theta),据此定义复数的 三角形式 (z = r(cos theta + isin theta))

利用 (cos)(sin) 的和角公式可得 (z_1z_2 = r_1r_2(cos(theta_1 + theta_2) + isin (theta_1 + theta_2)))。该等式称为 棣莫弗定理,它说明 复数相乘,模长相乘,辐角相加

  • 根据棣莫弗定理,我们得到对虚数单位 (i) 的一种直观理解:将一个复数 (z) 乘以 (i) 相当于将其 逆时针旋转 (frac {pi} 2) 弧度。实际上,考虑虚数单位 (i) 本身在复平面上的位置,发现就是 (1) 逆时针旋转 (frac {pi} 2) 度。一般地,有旋转的地方就有 (i) 存在,(i) 即旋转。推荐观看:虚数的来源 - Veritasium

1.3 单位根的定义

(r = 1) 时,(z = costheta + isintheta) 在单位圆上。此时根据棣莫弗公式有 (z ^ n = cos(ntheta) + sin(ntheta)),它在 复数旋转复数乘法 之间构建了桥梁:(z)(n) 次幂相当于从 ((1, 0)) 开始,以 (z) 的角度在单位圆上旋转 (n) 次得到的结果,称为将 (z) 旋转 (n) 次。

考虑将单位圆 (n) 等分(如下图),取任意 (n) 等分点 (P_k(0leq k ,将其旋转 (n) 次均得到 (1),因为跨过的 (frac 1 n) 单位圆弧数为 (n) 的倍数。这说明 (P_k ^ n = 1)

多项式 I:拉格朗日插值与快速傅里叶变换插图

探究 (P_k) 的表达式:因为它相当于从 (1) 开始在单位圆上旋转 (frac {2kpi} n) 弧度,因此 (P_k = cosleft(frac {2kpi} nright) + sinleft(frac {2kpi} nright))。我们称所有 (P_k)(n)单位根,将 (P_1) 记作 (omega_n),则 (P_k = P_1 ^ k = omega_n ^ k)

根据 (P_k ^ n = 1) 可知任意 (n) 次单位根 (omega_n ^ k) 均为 (x ^ n - 1 = 0) 的解。除特殊说明外,一般 (n) 次单位根直接代指 (omega_n),即从 (1) 开始逆时针方向的第一个单位根。

可以观看 3b1b 的视频 从 6:00 到 7:30 的部分加深对上述内容的理解。

  • 单位根的循环性:由 (omega_n ^ n = 1) 单位根的指数可对 (n) 取模。换言之,考虑从 (1) 开始不断乘以 (omega_n),我们将得到 (1, omega_n, omega_n ^ 2, cdots, omega_{n} ^ {n - 1}, omega_n ^ n = 1, omega_n ^ {n + 1} = omega_n, cdots),循环节为 (n)。想象从 (1) 开始不断旋转 (frac {2pi} n) 弧度,旋转 (n) 次后我们将落入循环。换言之,(omega_n ^ k = omega_n ^ {k + tn}(tin mathbb Z))
  • (omega_n ^ k = omega_{cn} ^ {ck}(c > 0))。对单位根有可视化认知后容易理解这一点。
  • (n) 为偶数时,将 (omega ^ k_n) 取反相当于将其逆时针(或顺时针)转半圈,所以 (-omega_n ^ k = omega_n ^ {kpm frac n 2}(2mid n))
  • 单位根的对称性:因为 (n) 次单位根将单位圆 (n) 等分,均匀分布在圆周,所以它们的重心就在原点,即 (sum_{i = 0} ^ {n - 1} omega_{n} ^ i = 0)
  • (gcd(k, n) = 1),则 (omega_n ^ k) 称为 本原单位根。所有本原单位根的 (0sim n - 1) 次幂互不相同。

1.4 单位根与原根

前置知识:原根,详见 初等数论学习笔记 I:同余相关

单位根和原根有极大的相似性,可以说原根就是模 (P) 下的整数域的单位根。

(n = varphi(P))(P) 存在原根 (g),则 (g ^ 0, g ^ 1, cdots, g ^ {n - 1}, g ^ n = g ^ 0, g ^ {n + 1} = g ^ 1) 这样的循环和 (n) 次单位根的循环一模一样。这使得在模 (P) 意义下涉及 (n) 次单位根运算时,可直接用原根 (g) 代替。进一步地,对于 (dmid n)(g ^ {frac n d}) 可直接代替模 (P) 意义下的 (d) 次单位根。

  • 单位根和原根都是对应域上一个大小为 (n = varphi(P))循环群生成元。它们均满足 (n) 次幂是对应域的单位元,且它们的 (0sim n - 1) 次幂互不相同。换言之,它们 同构

快速傅里叶变换 FFT 和快速数论变换 NTT 的联系恰在于此。

1.5 欧拉公式

前置知识:自然对数 (ln) 的底数 (e) 及其相关性质。

这部分为接下来可能用到的符号做一些补充。

欧拉公式指出

[e ^ {ix} = cos x + isin x
]

即单位圆上从 ((1, 0)) 开始旋转 (x) 弧度得到的复数,也即大小为 (x) 的角的终边与单位圆的交点。

欧拉公式的严谨证明超出了讨论范围,相关资料可以参考百度百科。我们给出一个对欧拉公式的感性理解方式,以加深读者对欧拉公式的直观印象与理解,来自 在 3.14 分钟内理解 (e ^ {ipi}) - 3Blue1Brown。这个视频相当有启发性。

根据 ((e ^ t)' = e ^ t),考虑 (e ^ t) 随着 (t) 增大在数轴上的取值。(e ^ t) 以每秒钟 (t) 均匀增加 (1) 的速度向数轴正方向移动,则 (e ^ t) 的速度就是它本身的位置。它的起始点为 (e ^ 0 = 1)

根据 ((e ^ {kt})' = ke ^ t),类似可知 (e ^ {kt}) 的速度等于 (k) 倍它本身的位置。

(k = i) 时,速度相当于将本身的位置逆时针旋转 (frac {pi} 2) 弧度,与位置垂直。这样,根据基础的高中物理知识,我们知道 (e ^ {it}) 随着 (t) 增大,在单位圆上做匀速圆周运动,且每秒移动 (1) 弧度。

这样,(e ^ {it}) 就等于模长为 (1),辐角为 (t) 的复数,即 (cos t + isin t)。这使得我们可以用 (r e ^ {itheta}) 描述模长为 (r),辐角为 (theta) 的复数,记号变得更简洁。

将该表示法应用至单位根,可知 (omega_n = e ^ {frac {2pi i} n})

读者应当在 (re ^ {itheta})(r(costheta + isin theta)) 及其在复平面上的位置建立直观联系,有助于理解接下来的内容。

带入 (t = pi),得到著名等式

[e ^ {ipi} = -1
]

带入 (t = 2pi = tau),得

[e ^ {itau} = 1
]

这说明对于任意 (kin mathbb Z)((e ^ {2pi i}) ^ {k + frac t n}) 相等恰对应 (omega_n ^ {kn + t}) 相等。

2. 多项式

2.1 基本概念

形如 (sum_{i = 0} ^ n a_ix ^ i)有限和式 称为 多项式,记作 (f(x) = sum_{i = 0} ^ n a_i x ^ i)。其中,(a_i) 称为 (i) 次项的 系数,也称 (x ^ i) 前的系数,记作 ([x ^ i]f(x))。超出最高次数 (n) 的系数 (a_i(i > n)) 视作 (0)

当项数无限时,和式形如 (sum_{i = 0} ^ {infty} a_ix ^ i),称为 形式幂级数,它暂时超出了我们的讨论范围。

  • 多项式 系数非零 的最高次项的次数称为该多项式的 ,也称次数,记作 (deg f)。如 (f(x) = sum_{i = 0} ^ n a_i x ^ i) 其中 (a_n neq 0),则 (f)(n) 次多项式,记作 (deg f = n)
  • 使得 (f(x) = 0) 的所有 (x) 称为多项式的
  • (a_i) 均为实数,则称 (f) 为实系数多项式。若 (a_i) 可以均为复数,则称 (f) 为复系数多项式。
  • 代数基本定理:任何非零一元 (n) 次复系数多项式恰有 (n) 个复数根。这些复数根可能重合。证明略。

2.2 系数表示法和点值表示法

(f(x) = sum_{i = 0} ^ n a_i x ^ i) 给出了所有 (i) 次项前的系数,这种描述多项式的方法称为 系数表示法

(x = x_i) 带入,得到 (y_i = f(x_i)),称 ((x_i, y_i))(f)(x_i) 处的点值。用若干点值 ((x_i, y_i)) 描述多项式的方法称为 点值表示法

考虑这两种表示法之间的联系。我们尝试探究在给定 (n) 个点值 ((x_0, y_0), (x_1, y_1), cdots, (x_{n - 1}, y_{n - 1})) 其中 (x_i) 互不相同时,所唯一确定的多项式的最高次数。

一个自然的猜测是 (n - 1),因为 (n - 1) 次多项式需要 (n) 个系数描述,具有 (n) 单位信息,而 (n) 个点值同样具有 (n) 单位信息。

证明即考虑直接带入,得到 (n) 元线性方程组:对于任意 (0leq j ,(sum_{i = 0} ^ {n - 1} a_ix_j ^ i = y_j)。该线性方程组的系数矩阵为

[begin{bmatrix}
1 & x_0 & x_0 ^ 1 & cdots & x_0 ^ {n - 1} \
1 & x_1 & x_1 ^ 1 & cdots & x_1 ^ {n - 1} \
1 & x_2 & x_2 ^ 1 & cdots & x_2 ^ {n - 1} \
vdots & vdots & vdots & ddots & vdots \
1 & x_{n - 1} & x_{n - 1} ^ 2 & cdots & x_{n - 1} ^ {n - 1} \
end{bmatrix}
]

(x_i) 互不相同,所以该范德蒙德矩阵的行列式非零,方程组有唯一解。类似的,从线性方程组的角度出发,容易证明 (n) 个点值不可能唯一确定 (n) 次或更高次的多项式。

综上,我们得到重要结论:(n) 个点值唯一确定的多项式最高次数为 (n - 1)

  • 从系数表示法转为点值表示法称为 求值(Evaluation)。
  • 从点值表示法转为系数表示法称为 插值(Interpolation)。

2.3 多项式运算

2.3.1 基本运算

(f(x) = sum_{i = 0} ^ n a_i x ^ i)(g(x) = sum_{i = 0} ^ m b_i x ^ i)

  • (h = f + g),则 (h(x) = (sum_{i = 0} ^ n a_i x ^ i) + (sum_{j = 0} ^ m b_j x ^ j) = sum_{i = 0} ^ {max(n, m)} (a_i + b_i) x ^ i),可知两多项式相加,对应系数相加,(deg (f + g) = max(deg f, deg g))
  • (h = f * g),则 (h(x) = (sum_{i = 0} ^ n a_i x ^ i)(sum_{j = 0} ^ m b_j x ^ j) = sum_{i = 0} ^ {n + m} (sum_{j = 0} ^ i a_jb_{i - j}) x ^ i),可知两多项式相乘,每两个系数相乘贡献至次数之和的系数,(deg(f * g) = deg f + deg g)

因此,在系数表示法下,计算两个多项式相加的复杂度为 (mathcal{O}(n + m)),计算两个多项式相乘的复杂度为 (mathcal{O}(nm))

  • 根据 ((f + g)(x) = f(x) + g(x)),可知两个多项式相加时,对应点值相加。
  • 根据 ((fg)(x) = f(x) g(x)),可知两个多项式相乘时,对应点值相乘。

因此,在点值表示法下,计算两个多项式相加需要 (max(n, m) + 1) 个点值,计算两个多项式相乘需要 (n + m + 1) 个点值,复杂度均为 (mathcal{O}(n + m))

  • (f * g)(fg) 表示多项式相乘,即进行加法卷积;用 (f cdot g) 表示多项式 点乘,即 对应系数相乘

在系数表示法下计算两个多项式相乘,我们首先将其转化为点值表示法,相乘后再转回系数表示法。时间复杂度 (mathcal{O}((n + m)log (n + m)))。相关算法将在第四章介绍。

3. 拉格朗日插值

在 2.2 小节我们得到结论:(n) 个点值唯一确定的多项式最高次数为 (n - 1)。那么,如何在点值表示法和系数表示法之间快速转换呢?

从系数表示法转为点值表示法,最常用的方法是直接带入,时间复杂度 (mathcal{O}(n ^ 2))(mathcal{O}(nlog ^ 2 n)) 的多项式多点求值则需要高级多项式技巧,此处不作介绍。

从点值表示法转为系数表示法,最直接的方法是高斯消元,时间复杂度 (mathcal{O}(n ^ 3))。接下来我们将介绍常用的拉格朗日插值法。

3.1 算法详解

拉格朗日插值的核心思想在于利用点值的可加性,每次只考虑一个点值,其它点值均视为 (0),由此构造 (n) 个多项式 (f_i(x)),使得它们在 (x_i) 处取值为 (y_i)(x_j(jneq i)) 处取值为 (0),则 (f = sum_{i = 0} ^ {n - 1} f_i)中国剩余定理 使用了类似的思想。

为了让 (f_i(x_j) = 0 (ineq j))(f_i) 应包含 (x - x_j) 作为因式,因此设 (f_i(x) = prod_{i ne j} (x - x_j))。但是此时 (f_i(x_i)) 不一定等于 (y_i),我们发现可以调整 (f_i) 的系数,给 (f_i) 乘上 (frac {y_i} {f_i(x_i)})。综上,我们得到 (f_i) 的表达式

[f_i(x) = y_i prod_{j neq i} frac {x - x_j} {x_i - x_j}
]

(f_i) 求和得 (f),我们得到拉格朗日插值公式

[f(x) = sumlimits_{i = 0} ^ {n - 1} y_i prodlimits_{j neq i} frac {x - x_j} {x_i - x_j}
]

为得到 (f) 的各项系数,只需 (mathcal{O}(n ^ 2)) 求出 (F(x) = prod_{i = 0} ^ {n - 1} (x - x_i)),对每个 (i) 暴力 (mathcal{O}(n))(F) 除掉一次式 (x - x_i) 算出 (frac {F(x)} {x - x_i}) 的各项系数,再乘以 (frac {y_i} {prod_{j neq i} x_i - x_j}) 得到 (f_i),则 (f = sum_{i = 0} ^ {n - 1} f_i)。时间复杂度 (mathcal{O}(n ^ 2))

通常情况下,题目要求我们求出 (F(x)) 在给定某个 (x) 处的取值,此时我们不把 (x) 看做函数的元,而是直接将 (x) 带入上式。时间复杂度仍为 (mathcal{O}(n ^ 2))

多项式快速插值在 (mathcal{O}(nlog ^ 2 n)) 的时间内将点值表示法转化为系数表示法,这超出了我们的讨论范围。

3.2 连续取值插值

当给定点值横坐标为连续整数时,我们有快速单点插值的方法。

(0, 1, cdots, n - 1)(x_i = i) 为例:

[begin{aligned}
f(x) = sumlimits_{i = 0} ^ {n - 1} y_i prodlimits_{j neq i} frac {x - j} {i - j} \
end{aligned}
]

分子是 (prod (x - i)) 挖掉一个项,经典维护前缀后缀积。设 (p_i = prod_{j = 0} ^ {i - 1} x - j)(s_i = prod_{j = i + 1} ^ {n - 1} x - j)

分母对于 (i > j)(prod_{j = 0} ^ {i - 1} (i - j) = i!)。对于 (i ,(prod_{j = i + 1} ^ {n - 1} (i - j) = (-1)(-2) cdots (i - n + 1)),将所有负号提出来,得 ((-1) ^ {n - i + 1}(i - n + 1)!)

因此

[f(x) = sumlimits_{i = 0} ^ {n - 1} y_i frac {p_is_i} {(-1) ^ {n - i + 1} i! (i - n + 1)!}
]

预处理阶乘逆元,时间复杂度 (mathcal{O}(n))

3.3 应用

  • 计算范德蒙德方阵的逆矩阵,详见 4.4.1 小节。

3.4 例题

P4781 【模板】拉格朗日插值

#include 
using namespace std;
constexpr int N = 2e3 + 5;
constexpr int mod = 998244353;
int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % mod;
    a = 1ll * a * a % mod, b >>= 1;
  }
  return s;
}
int n, k, ans, x[N], y[N];
int main() {
  cin >> n >> k;
  for(int i = 1; i > x[i] >> y[i];
  for(int i = 1; i 

4. 快速傅里叶变换

快速傅里叶变换(Fast Fourier Transform,FFT)是多项式算法的根基。想要真正理解它,必须先真正理解单位根,还需要对线性代数有基本了解。

推荐观看:

4.1 求值的本质

(f(x) = sum_{i = 0} ^ {n - 1} a_i x ^ i),将 (x_0) 带入,得 (f(x_0) = sum_{i = 0} ^ {n - 1} a_i x_0 ^ i)。考虑将其写成向量内积(点积)的形式:

[f(x_0) = sum_{i = 0} ^ {n - 1} a_i x_0 ^ i =
begin{bmatrix}
x_0 ^ 0 & x_0 ^ 1 & cdots & x_0 ^ {n - 1}
end{bmatrix}
times
begin{bmatrix}
a_0 \ a_1 \ vdots \ a_{n - 1}
end{bmatrix}
]

这样,如果有 (x_0, x_1, cdots, x_{m - 1}) 需要求值,整个过程就可以写成 (mtimes n) 维矩阵乘以 (n) 维列向量的形式:

[begin{bmatrix}
x_0 ^ 0 & x_0 ^ 1 & cdots & x_0 ^ {n - 1} \
x_1 ^ 0 & x_1 ^ 1 & cdots & x_1 ^ {n - 1} \
vdots & vdots & ddots & vdots \
x_{m - 1} ^ 0 & x_{m - 1} ^ 1 & cdots & x_{m - 1} ^ {n - 1} \
end{bmatrix}
times
begin{bmatrix}
a_0 \ a_1 \ vdots \ a_{n - 1}
end{bmatrix}
=
begin{bmatrix}
f(x_0) \ f(x_1) \ vdots \ f(x_{m - 1})
end{bmatrix}
]

左侧矩阵就是著名的 范德蒙德矩阵

(m = n) 时为范德蒙德方阵,(x_i) 互不相同时其逆存在,帮助我们快速从点值表示法转回系数表示法。(m > n) 时任取 (x_i) 互不相同的 (n + 1) 行可以求逆。(m 时无法还原系数。这体现出 (n + 1) 个点值唯一确定最高次数不超过 (n) 的多项式。

朴素计算求值的复杂度为 (mathcal{O}(nm)),因为带入求值一次的复杂度为 (mathcal{O}(n))。快速傅里叶变换即在离散傅里叶变换基础上通过选取合适的 (x_i),使得可以快速求值。

4.2 离散傅里叶变换

在介绍 FFT 之前,我们先给出离散傅里叶变换(Discrete Fourier Transform,DFT)的概念。

DFT 在工程中是将离散信号从时域转为频域的过程。碰巧的是,其表达式刚好可以用来对多项式进行多点求值,只不过这些点值是固定的 单位根 处的点值,但对于求值做多项式乘法已经足够了。

DFT 将一个长度为 (N) 的复数序列 (x_0, x_1, cdots, x_{N - 1}) 通过如下公式转化为另一个长为 (N) 的复数序列 (X_0, X_1, cdots, X_{N - 1})

[X_k = sum_{n = 0} ^ {N - 1} x_n e ^ {-frac {2pi i} Nkn}
]

也即

[X_k = sum_{n = 0} ^ {N - 1} x_n omega_N ^ {-kn}
]

设多项式 (f(x) = sum_{n = 0} ^ {N - 1} x_n x ^ i),易知

[X_k = sum_{n = 0} ^ {N - 1} x_n(omega_N ^ {-k}) ^ n = f(omega_N ^ {-k})
]

根据上式,离散傅里叶变换可以看成多项式 (f(x) = sum_{n = 0} ^ {N - 1} x_nx ^ i)(N) 个单位根处求值。

没看懂?没关系。接下来我们将从另一角度出发,得到一个差别微小但更容易理解的算法 —— FFT。

4.3 算法详解

首先我们得搞清楚,DFT 是一个变换,而 FFT 是用于实现 DFT 的算法。在 FFT 的推导过程中,其所实现的变换和 DFT 变换有细微的差别,因此笔者也用 FFT 表示该算法实现的变换。

按理说学习一个算法时需要搞清楚它的用途,但如果直接从 DFT 角度入手尝试简化流程,那么为了让动机看起来自然,我们还要先学习 DFT 的实际含义,这涉及到大量前置知识。

但是,从计算机科学界尤为重要且为各位 OIer 熟知的多项式乘法的优化出发,我们通过自然推理得到一个算法,而该算法恰好可以快速实现 DFT(有一些细小的差别,详见 4.3.5)。

我们明确接下来的目标:给定次数 (n - 1) 的多项式 (f(x) = sum_{i = 0} ^ {n - 1} a_i x ^ i(a_{n - 1} neq 0)),快速求出它的至少 (n) 个点值。

下文中,(f(x)) 也被视为关于 (x)(n - 1) 次函数。

4.3.1 简化情况

解决一个普遍性的问题,最重要的思想就是 从简化情况入手,分析问题的性质

函数的性质无非就那几个:奇偶性,单调性,周期性等。一般函数没有单调性和周期性,但总可以表示为一个偶函数和一个奇函数之和,这一定是突破点。

  • 对于偶函数 (f(x)),即所有奇数次项系数为 (0),带入两个相反数 (w)(-w) 时,它们的值相等。
  • 对于奇函数 (f(x)),即所有偶数次项系数为 (0),带入两个相反数 (w)(-w) 时,它们的值互为相反数。

因此,将任意多项式 (f(x)) 拆成偶函数 (f_e(x)) 和奇函数 (f_o(x)) 之和,则

[begin{cases}
f(x) = f_e(x) + f_o(x) \
f(-x) = f_e(x) - f_o(x)
end{cases}
]

选择 (lceilfrac n 2rceil) 对两两互为相反数的值 ((x_i, -x_i)) ,求出所有 (x_i)(f_e(x))(f_o(x)) 处的取值。

不妨设 (n) 为偶数,则 (f_e)(n - 2) 次多项式,(f_e)(n - 1) 次多项式,本质上依然相当于求出 (n - 1) 次多项式的 (n) 个点值,对时间复杂度没有优化。

但是 (f_e)(f_o) 的项数减半,尝试利用该性质。

因为 (f_e) 只有偶数次项 (a_0x ^ 0 + a_2x ^ 2 + cdots),故考虑换元 (u = x ^ 2),则 (f_e(u) = a_0u ^ 0 + a_2 u ^ 1 + cdots)。换言之,我们设 (f'_e(x) = a_0x ^ 0 + a_2x ^ 1 + cdots),则 (f_e(x) = f'_e(x ^ 2))

同理,从 (f_o) 中提出一个 (x),设 (f'_o(x) = a_1x ^ 0 + a_3x ^ 1 + cdots),则 (f_o(x) = xf'_o(x ^ 2))

因此,

[begin{cases}
f(x) = f'_e(x ^ 2) + xf'_o(x ^ 2) \
f(-x) = f'_e(x ^ 2) - xf'_o(x ^ 2)
end{cases}
]

这样才是真正意义上的规模减半。若问题可递归,则时间复杂度为 (T(n) = 2Tleft(frac n 2right) + mathcal{O}(n)),解得 (T(n) = mathcal{O}(nlog n))

4.3.2 引入单位根

问题来了,如何保证递归得到的问题也满足两两互为相反数呢?

考虑一开始的 ((x_i, -x_i)),这说明存在 (i') 使得 (x_i ^ 2 = -x_{i'} ^ 2),它们互不相同但它们的 (4) 次方相等。

进一步推演,因为问题会递归 (w = lceillog_2 nrceil) 层,所以我们需要找到 (k = 2 ^ w) 个互不相等的 (x),但它们的 (k) 次幂相等。这个相等的 (k) 次幂可以任意选择,方便起见设为 (1),则 (x ^ k = 1)(x) 为所有 (k) 次单位根。

逆推得到结果后,我们再顺着检查一遍:初始令 (x) 为所有 (k) 次单位根,每层递归会将这些数平方,得到所有 (frac k 2, frac k 4 cdots) 次单位根。因为 (frac k 2, frac k 4, cdots) 都是偶数,所以它们可以两两正负配对,直到递归 (w) 层的平凡情况:(frac k {2 ^ w} = 1) 次单位根即 (x = 1)

由此可得通常使用(补齐到 (2) 的幂)的快速傅里叶变换的范德蒙德方阵形式:

[begin{bmatrix}
(omega_n ^ 0) ^ 0 & (omega_n ^ 0) ^ 1 & cdots & (omega_n ^ 0) ^ {n - 1} \
(omega_n ^ 1) ^ 0 & (omega_n ^ 1) ^ 1 & cdots & (omega_n ^ 1) ^ {n - 1} \
vdots & vdots & ddots & vdots \
(omega_n ^ {n - 1}) ^ 0 & (omega_n ^ {n - 1}) ^ 1 & cdots & (omega_n ^ {n - 1}) ^ {n - 1} \
end{bmatrix}
times
begin{bmatrix}
a_0 \ a_1 \ vdots \ a_{n - 1}
end{bmatrix}
=
begin{bmatrix}
f(omega_n ^ 0) \ f(omega_n ^ 1) \ vdots \ f(omega_n ^ {n - 1})
end{bmatrix}
]

4.3.3 递归公式

得到大致框架后,我们具体地描述整个算法流程:

首先将 (n) 补齐到不小于 (n) 的最小的 (2) 的幂,即 (2 ^ {lceil log_2 nrceil})

设当前需要求值的多项式 (f) 长度为 (L(L = 2 ^ w, win mathbb N)),若 (L = 1) 则直接返回 (a_0)。否则我们需求出 (f(x)) 在所有 (L) 次单位根 (omega_L ^ k(0leq k 处的取值。

(f(x)) 写成 (f_e(x ^ 2) + xf_o(x ^ 2)),则对于 (0leq k ,有

[begin{aligned}
f(omega_L ^ k) & = f_e(omega_L ^ {2k}) + omega_L ^ k f_o(omega_L ^ {2k}) \
f(omega_L ^ {k + frac L 2}) & = f_e(omega_L ^ {2k + L}) + omega_L ^ {k + frac L 2} f_o(omega_L ^ {2k + L})
end{aligned}
]

根据单位根的性质(在单位根部分有介绍):

  • (omega_n ^ k = omega_{2n} ^ {2k})
  • (omega_n ^ k = omega_n ^ {k + tn}(tin mathbb Z))
  • (omega_{n} ^ k = -omega_{n} ^ {k + frac n 2} (2mid n))

[begin{aligned}
f(omega_L ^ k) & = f_e(omega_{frac L2} ^ {k}) + omega_L ^ k f_o(omega_{frac L 2} ^ {k}) \
f(omega_L ^ {k + frac L 2}) & = f_e(omega_{frac L2} ^ {k}) - omega_L ^ k f_o(omega_{frac L 2} ^ {k})
end{aligned}
]

这样,我们只需求出 (frac L 2) 次多项式 (f_e)(f_o) 在所有 (frac L 2) 次单位根处的取值。

4.3.4 蝴蝶变换

递归处理比较慢,我们希望像位运算卷积一样通过递推实现整个过程。

因为在边界条件下,每个位置的取值与其对应的系数相关,所以递归转递推的关键在于考察系数的变化。

考虑 (n = 8) 的情况,初始为 ((a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7))

进行第一层递归时,将 (f_e) 的系数写在左半部分,(f_o) 的系数写在右半部分,得 ((a_0, a_2, a_4, a_6), (a_1, a_3, a_5, a_7))

进行第二层递归时,类似地将每个子问题的 (f_e)(f_o) 的系数分别写在左右两侧,得 ((a_0, a_4), (a_2, a_6), (a_1, a_5), (a_3, a_7))

进行第三层递归时,共有 (4) 个大小为 (2) 的子问题,且进行上述操作后系数的位置不发生改变。

我们看到,如果一个系数在规模为 (2n) 的问题中的位置为 (p),若 (p) 是偶数,则递归至左半部分,若 (p) 是奇数,则递归至右半部分,且在规模为 (n) 的问题中的位置为 (leftlfloor frac p 2rightrfloor)

进一步地,一个系数在第 (i) 次递归时的方向决定了它最终下标在二进制下第 (w - i(n = 2 ^ w)) 位的取值。若向左递归则为 (0),向右递归则为 (1)。而它递归的方向又受到 (leftlfloor frac p {2 ^ {i - 1}}rightrfloor) 的奇偶性的影响,即 (p) 在二进制下第 (i) 位的取值,若为 (0) 则向左递归,为 (1) 则向右递归。

这就说明,(p) 在二进制下第 (i) 位的取值,等于它对应的系数的最终下标在二进制下第 (w - i) 位的取值。

因此我们断言,系数 (a_p)(w) 次递归后的下标等于 (p) 二进制翻转后的值,设为 (r_p)。这里翻转指翻转第 (0sim w - 1) 位的值。

(r_p) 可递推求得:(r_0 = 0)。对于 (r_i(i > 0)),先右移一位得到 (i' = leftlfloor frac i 2rightrfloor),则 (r_i) 的低 (w - 1) 位(第 (0sim w - 2) 位)的值可由 (r_{i'}) 右移一位得到。第 (w - 1) 位的值只需检查 (i) 的奇偶性。因此,(r_i = leftlfloor frac {r_{i'}}{2} rightrfloor + frac n 2(ibmod 2)),其中 (i' = lfloor frac i 2rfloor)

因此,在算法一开始先对系数数组 (a) 执行 蝴蝶变换,即同时令 (a_i to a_{r_i}),然后类似 FWT,枚举问题规模 (2d(1le d ,枚举每个子问题 (2di(0leq 2di ,再枚举当前子问题的所有对应位置 ((x = 2di + k, y = 2di + k + d)(0leq k 执行变换,即设 (x)(y) 对应位置的当前值为 (f_x)(f_y),则 (f'_x = f_x + omega_{2d} ^ k f_y)(f'_y = f_x - omega_{2d} ^ k f_y)

进一步地,根据 (r) 的定义,我们有 (r_{r_i} = i)。因此,执行蝴蝶变换只需对所有 (i 执行 (mathrm{swap}(a_i, a_{r_i}))

这就是 FFT,整个过程称为 对多项式 (f) 做长度为 (n) 的快速傅里叶变换,时间复杂度 (mathcal{O}(nlog n))。代码在 4.5 小节。

4.3.5 DFT 和 FFT

对比 DFT 和 FFT 矩阵:

[mathcal {F_D} =
begin{bmatrix}
(omega_N ^ 0) ^ 0 & (omega_N ^ 0) ^ 1 & cdots & (omega_N ^ 0) ^ {N - 1} \
(omega_N ^ {-1}) ^ 0 & (omega_N ^ {-1}) ^ 1 & cdots & (omega_N ^ {-1}) ^ {N - 1} \
vdots & vdots & ddots & vdots \
(omega_n ^ {-(N - 1)}) ^ 0 & (omega_n ^ {-(N - 1)}) ^ 1 & cdots & (omega_N ^ {-(N - 1)}) ^ {N - 1} \
end{bmatrix}
\
mathcal {F_F} =
begin{bmatrix}
(omega_n ^ 0) ^ 0 & (omega_n ^ 0) ^ 1 & cdots & (omega_n ^ 0) ^ {n - 1} \
(omega_n ^ 1) ^ 0 & (omega_n ^ 1) ^ 1 & cdots & (omega_n ^ 1) ^ {n - 1} \
vdots & vdots & ddots & vdots \
(omega_n ^ {n - 1}) ^ 0 & (omega_n ^ {n - 1}) ^ 1 & cdots & (omega_n ^ {n - 1}) ^ {n - 1} \
end{bmatrix}
]

可以发现 DFT 和 FFT 基本一致,它们的差别在于:

  • 朴素 FFT 要求 (n)(2) 的幂,但 DFT 序列长度可以是任意正整数。
  • DFT 和 FFT 带入单位根的顺序是相反的。

注意这些细微差别,不要把 DFT 和 FFT 搞混了。

4.3.6 循环卷积

4.4 离散傅里叶逆变换

离散傅里叶逆变换(Inverse Discrete Fourier Transform,IDFT)可以视为单位根处插值的过程。即给出 (n = 2 ^ w) 个在所有 (n) 次单位根处的点值 (P_k = (omega_n ^ k, f(omega_n ^ k))(0leq k ,要求还原 (f) 的各项系数,其中 (f) 的次数不大于 (n - 1)

类似地,IDFT 和 IFFT 之间也存在一些差异。想必各位读者根据之前的内容已经可以猜出这种差异是什么了。

4.4.1 范德蒙德方阵逆

考虑范德蒙德方阵

[mathcal A = begin{bmatrix}
x_0 ^ 0 & x_0 ^ 1 & cdots & x_0 ^ {n - 1} \
x_1 ^ 0 & x_1 ^ 1 & cdots & x_1 ^ {n - 1} \
vdots & vdots & ddots & vdots \
x_{n - 1} ^ 0 & x_{n - 1} ^ 1 & cdots & x_{n - 1} ^ {n - 1} \
end{bmatrix}
]

这玩意怎么求逆?我们考虑最暴力的方法:拉格朗日插值!

因为范德蒙德方阵可以看成多项式多点求值,根据

[begin{bmatrix}
x_0 ^ 0 & x_0 ^ 1 & cdots & x_0 ^ {n - 1} \
x_1 ^ 0 & x_1 ^ 1 & cdots & x_1 ^ {n - 1} \
vdots & vdots & ddots & vdots \
x_{n - 1} ^ 0 & x_{n - 1} ^ 1 & cdots & x_{n - 1} ^ {n - 1} \
end{bmatrix}
times
begin{bmatrix}
a_0 \ a_1 \ vdots \ a_{n - 1}
end{bmatrix}
=
begin{bmatrix}
f(x_0) \ f(x_1) \ vdots \ f(x_{n - 1})
end{bmatrix}
]

再结合拉格朗日插值公式

[f(x) = sumlimits_{i = 0} ^ {n - 1} f(x_i) prodlimits_{j neq i} frac {x - x_j} {x_i - x_j}
]

可知从 (f(x_j)) 贡献到 (a_i) 的系数为 ((mathcal{A} ^ {-1})_{ij} = [x ^ i] prod_{kneq i} frac {x - x_k} {x_j - x_k})

这就是范德蒙德方阵逆当中每个元素的表达式。

4.4.2 算法介绍

很显然,(f) 在经过快速傅里叶变换后,再进行快速傅里叶逆变换,仍得到 (f)

因此,只需对快速傅里叶变换的矩阵求逆,即可得到快速傅里叶逆变换的矩阵。

(omega) 表示 (omega_n),则

[mathcal F =
begin{bmatrix}
(omega ^ 0) ^ 0 & (omega ^ 0) ^ 1 & cdots & (omega ^ 0) ^ {n - 1} \
(omega ^ 1) ^ 0 & (omega ^ 1) ^ 1 & cdots & (omega ^ 1) ^ {n - 1} \
vdots & vdots & ddots & vdots \
(omega ^ {n - 1}) ^ 0 & (omega ^ {n - 1}) ^ 1 & cdots & (omega ^ {n - 1}) ^ {n - 1} \
end{bmatrix}
]

((mathcal F ^ {-1})_{ij} = [x ^ i] prod_{kneq j} frac {x - omega ^ k} {omega ^ j - omega ^ k})

分子和分母均形如 (g(x) = frac {prod_{0leq k ,我们研究该函数的性质。

首先,因为 (omega ^ k(0leq k 为 (x ^ n - 1 = 0)(n) 个互不相同的根,根据因式定理,(prod_{0leq k 就等于 (x ^ n - 1)。感性理解:将 (prod_{0leq k 展开,再应用单位根的 对称性

模拟短除法 (frac {x ^ n - 1} {x - omega ^ j}),可知

[g(x) = x ^ {n - 1} + omega ^ jx ^ {n - 2} + (omega ^ j) ^ 2x ^ {n - 3} + cdots + (omega ^ j) ^ {n - 1} x ^ 0
]

因此

[(mathcal F ^ {-1})_{ij} = frac {[x ^ i] g(x)} {g(omega ^ j)} = frac {(omega ^ j) ^ {n - 1 - i}} {n(omega ^ j) ^ {n - 1}} = frac {(omega ^ {-j}) ^ {i}omega ^ {-j}} {nomega ^ {-j}} = frac {omega ^ {-ij}} {n}
]

这就说明

[mathcal F ^ {-1} =
frac 1 n
begin{bmatrix}
(omega ^ {-0}) ^ 0 & (omega ^ {-0}) ^ 1 & cdots & (omega ^ {-0}) ^ {n - 1} \
(omega ^ {-1}) ^ 0 & (omega ^ {-1}) ^ 1 & cdots & (omega ^ {-1}) ^ {n - 1} \
vdots & vdots & ddots & vdots \
(omega ^ {-(n - 1)}) ^ 0 & (omega ^ {-(n - 1)}) ^ 1 & cdots & (omega ^ {-(n - 1)}) ^ {n - 1} \
end{bmatrix}
]

因此,对一个序列做 IFFT,只需将 FFT 递归公式里面的 (omega_L ^ k) 换成 (omega_L ^ {-k}),并在最后将所有数除以 (n) 即可。

这就引出了 IDFT 公式及其对应矩阵:

[x_n = frac 1 N sum_{k = 0} ^ {N - 1} X_k e ^ {frac {2pi i} Nkn} = sum_{k = 0} ^ {N - 1} X_k omega_N ^ {kn} \
mathcal {F_D} ^ {-1} =
frac 1 N
begin{bmatrix}
(omega_N ^ 0) ^ 0 & (omega_N ^ 0) ^ 1 & cdots & (omega_N ^ 0) ^ {N - 1} \
(omega_N ^ 1) ^ 0 & (omega_N ^ 1) ^ 1 & cdots & (omega_N ^ 1) ^ {N - 1} \
vdots & vdots & ddots & vdots \
(omega_N ^ {N - 1}) ^ 0 & (omega_N ^ {N - 1}) ^ 1 & cdots & (omega_N ^ {N - 1}) ^ {N - 1} \
end{bmatrix}
]

4.5 快速多项式乘法

通过 FFT 和 IFFT,我们可以在 (mathcal{O}(nlog n)) 的时间内实现 (n - 1) 次多项式在系数表示法和点值表示法之间的转换。

计算两个次数分别为 (n - 1)(m - 1) 的多项式 (a, b) 相乘,设结果为 (c),则 (c)(n + m - 2) 次多项式,我们需要 (n + m - 1) 个点值才能确定它。根据点值的性质,首先找到不小于 (n + m - 1) 的最小的 (2) 的幂 (L),对系数表示法的 (a, b) 分别做长度为 (L) 的 FFT,将对应点值相乘得到 (hat c),再对 (hat c) 做 IFFT 还原 (c)

多项式 I:拉格朗日插值与快速傅里叶变换插图1


首先我们需要实现一个复数类,它支持复数的加减乘运算。也可以使用 C++ 自带复数类 std::complex,用法见 CPP reference

FFT 和 IFFT 大体上一致,只有一些细小差别。我们可以类似实现 FWT 和 IFWT 那样,用同一个函数实现它们,并用一个参数区别。

等式 (omega_n = cos(frac {2pi} {n}) + isin(frac {2pi} {n})) 帮助我们求出 (n) 次单位根。

  • 注意浮点数的精度。当多项式系数的绝对值较大时,需使用 long double 甚至 5.2 小节的任意模数卷积。

模板题 P3803 代码:

#include 
using namespace std;

constexpr int N = 1 > 1] >> 1) + (i & 1 ? L >> 1 : 0);
    if(i > n >> m;
  for(int i = 0; i > f[i].a;
  for(int i = 0; i > g[i].a;
  int L = 1;
  while(L 

4.6 快速数论变换

前置知识:原根和阶。

我们将 DFT 的数值范围从复数域推广至任意存在 (n) 次单位根 (alpha) 的环 (R),其中 (alpha) 满足 (alpha ^ k(0leq k 互不相同,对 (x_0, x_1, cdots, x_{n - 1}) DFT 得 (X_0, X_1, cdots, X_{n - 1}),则

[X_i = sum_{j = 0} ^ {n - 1} x_j alpha ^ {ij}
]

也可以视作 (X_i = f(alpha ^ i)),其中 (f(t) = sum_{j = 0} ^ {n - 1} x_j t ^ j)

类似可知 IDFT

[x_j = frac 1 n sum_{i = 0} ^ {n - 1} X_i alpha ^ {-ij}
]

即 DFT 和 IDFT 对序列进行的变换的本质抽象。

4.6.1 算法介绍

快速数论变换即在模 (p) 意义下的整数域 (F = mathbb Z / p) 进行的快速傅里叶变换。

设变换长度为 (n),则需存在 (n) 次单位根 (a) 满足 (delta_p(a) = n)。大部分题目的模数 (p) 满足:

  • (p) 为质数。
  • (2 ^ kmid p - 1),其中 (2 ^ k) 不小于最大的可能的 NTT 长度。

第一条性质保证 (p) 存在原根 (g)(n) 可逆,第二条性质保证存在 (n) 次单位根。注意它不是充要条件,只是充分条件。

根据原根的性质,(delta_p(g) = p - 1),即 (g)(0sim p - 2) 次幂互不相同,则 (g)(F) 上的性质和复数域上 (p - 1) 次单位根的性质完全一致:(g ^ k(0leq k 和 (omega_{p - 1} ^ k(0leq k 形成的域是同构的,(g ^ k) 等价于 (omega_{p - 1} ^ k)

因此,(q = g ^ {frac {p - 1} {n}}) 等价于 (omega_{p - 1} ^ {frac {p - 1} n})(omega_n),它的 (0sim n - 1) 次幂互不相同,即 (delta_p(q) = n),也可以通过阶的性质 (delta_p(g ^ k) = frac {delta_p(g)} {gcd(delta_p(g), k)}) 说明。

常见 NTT 模数有:

  • (998244353 = 7times 17 times 2 ^ {23} + 1),有原根 (3)。它可以用来做长度不超过 (2 ^ {23}) 的 NTT,也是最常见的模数。
  • (1004535809 = 479 times 2 ^ {21} + 1),有原根 (3)
  • (469762049 = 7 times 2 ^ {26} + 1),有原根 (3)

如果不是常见模数,我们需要检验 (p) 是否为形如 (q2 ^ k + 1) 的质数,(2 ^ k) 是否足够大,再运用原根相关的知识枚举并判定找到任意一个原根,就可以做 NTT 了。

注意以上只是模数 (p) 可 NTT 的充分条件,它的更弱的条件是存在 (delta_{p}(a) = n)(n ^ {-1})。如果 NTT 是正解的一部分,那么一个合格的出题人应将 (p) 设为常见模数,因为模数不是考察重点。

4.6.2 代码实现

朴素 NTT 已经比 FFT 快了不少,因为复数运算很耗时。

接下来加入一些常数优化:

  • 预处理 (2d) 次单位根的 (0sim d - 1) 次幂,这样对每个 (i) 枚举 (j) 的时候就不需要重复计算单位根的幂。
  • unsigned long long(16) 次模一次的技巧,减少取模次数。类似技巧也用于优化矩阵乘法。

综上,写出如下代码(依然是 模板题):尽管题目没有要求取模,但可视为在模 (p) 意义下进行多项式乘法,其中 (p) 大于答案的每一项。

#include 
using namespace std;

using ull = unsigned long long;
constexpr int N = 1 >= 1;
  }
  return s;
}

int n, m, r[N], f[N], g[N], h[N];
void FFT(int L, int *a, bool type) {
  static ull f[N], w[N];
  for(int i = 0; i > 1] >> 1) | (i & 1 ? L >> 1 : 0)];
  for(int d = 1; d > n >> m;
  for(int i = 0; i > f[i];
  for(int i = 0; i > g[i];
  int L = 1;
  while(L 

5. 应用与技巧

FFT 和 NTT 是所有快速多项式操作的基础。

设多项式 (f) DFT 得到 (hat f),也记为 (operatorname {DFT}(f))。可以发现,DFT 是线性变换,因此它具有 线性性,这是它最重要且最常用的一个性质:

  • (c operatorname {DFT}(f) + d operatorname {DFT}(g) = operatorname {DFT} (cf + dg))

5.1 常数优化

在分析变换次数的时候,一般不区分 DFT 和 IDFT。

5.1.1 三次变两次优化

计算两 实系数 多项式 (f, g) 相乘。

根据 ((a + bi) ^ 2 = (a ^ 2 - b ^ 2) + 2abi),将 (f + gi) 平方后取出虚部再除以 (2) 即可。

这样,三次 DFT 变成了两次 DFT,对常数有显著优化。代码

这个优化被接下来稍复杂的技巧完全包含,它们的核心思想都是利用实系数的性质。

5.1.2 合并两次实系数 DFT

同时计算两 实系数 多项式 (f, g) 的 DFT。

类似地,我们设 (u = f + ig)(v = f - ig)。先求出 (u) 的 DFT (hat u),考虑能否利用 (hat u) 直接求出 (hat v)

在给出做法之前,我们还要引入一些复数相关的概念:

  • 定义:对于两个复数 (z_1)(z_2),若它们实部相等,虚部互为相反数,则称 (z_1, z_2) 互为 共轭复数(z_2)(z_1) 的共轭,(z_1)(z_2) 的共轭。
  • 表示:复数 (z) 的共轭记作 (overline {z})(operatorname {conj}(z))
  • 运算性质:两个共轭复数的辐角互为相反数,即 (arg z_1 = -arg z_2)。根据棣莫弗定理,可知 复数积的共轭,等于它们共轭的积。这样理解:求共轭相当于把整个复平面沿 (x) 轴翻转,求积再翻转等价于翻转再求积。
  • 易知复数和的共轭等于它们共轭的和,且一个复数的共轭的共轭等于它本身。

首先,(v(omega_n ^ k) = sum_{i = 0} ^ {n - 1} v_i(omega_n ^ k) ^ i)。为了让它和 (u) 产生联系,因为 (f, g) 为实系数多项式,所以 (u)(v) 的各项系数互为共轭,我们对整个结果求两次共轭,并将一次共轭根据共轭的性质下放至 (v_i)(omega_n ^ k)

[begin{aligned}
hat v_k
& = v(omega_n ^ k) \
& = sum_{i = 0} ^ {n - 1} v_i(omega_n ^ k) ^ i \
& = operatorname {conj} left( operatorname {conj} left(sum_{i = 0} ^ {n - 1} v_i(omega_n ^ k) ^ i right) right) \
& = operatorname {conj} left( sum_{i = 0} ^ {n - 1} operatorname {conj}(v_i(omega_n ^ k) ^ i) right) \
& = operatorname {conj} left( sum_{i = 0} ^ {n - 1} operatorname {conj}(v_i)operatorname {conj}(omega_n ^ k) ^ i right) \
& = operatorname {conj} left( sum_{i = 0} ^ {n - 1} operatorname {conj}(v_i)operatorname {conj}(omega_n ^ k) ^ i right) \
& = operatorname {conj} left( sum_{i = 0} ^ {n - 1} u_i (omega_n ^ {n - k}) ^ i right) \
& =
begin{cases}
operatorname {conj} (hat u ^ {n - k}) & (1leq k

求得 (hat v) 之后,因为 (hat u = hat f + ihat g)(hat v = hat f - ihat g),所以 (hat f = frac {hat u + hat v} {2})(hat g = frac {hat u - hat v} {2i} = frac {(hat v - hat u)i} {2})

5.1.3 合并两次实系数 IDFT

这里实系数指 IDFT 后的各项系数为实数。如果是 IDFT 前的各项系数为实数,直接使用上一小节的技巧即可。

设需要 IDFT 的两个多项式为 (hat f)(hat g)。计算 (operatorname {IDFT}(hat f)),各项系数均为实数,虚部信息被浪费了。考虑计算 (operatorname {IDFT}(hat f + ihat g)),则各项系数的实部即 (f) 的系数,虚部即 (g) 的系数。

5.2 任意模数卷积

给定两多项式 (f, g),在系数模 (p) 意义下求出它们的卷积 (h = f * g)

(p) 不是 NTT 模数,我们不能朴素地使用 FFT 求解该问题,因为 (h) 的系数可达 (n p ^ 2)。取 (n = 10 ^ 6)(p = 10 ^ 9),则 (np ^ 2 = 10 ^ {24})long double 也无法接受。

接下来介绍两种常见的解决该问题的方法。它们也被称为 MTT(非正式),得名于毛啸 2016 年的集训队论文。

5.2.1 拆系数 FFT

(B = sqrt p),将所有系数表示为 (kB + r(0leq r 的形式,得到四个多项式 (f = f_0B + f_1)(g = g_0B + g_1),计算它们两两相乘的结果,则 (fg = (f_0g_0) B ^ 2 + (f_0g_1 + f_1g_0)B + f_1g_1)

显然有一个四次 DFT 和三次 IDFT 的朴素做法:对 (f_0, f_1, g_0, g_1) 进行 DFT,(hat f_0cdot hat g_0, hat f_0cdot hat g_1 + hat f_1cdot hat g_0, hat f_1 cdot hat g_1) 进行 IDFT。

使用 6.1.2 和 6.1.3 的技巧,可以做到两次 DFT 和两次 IDFT。系数值域 (nB ^ 2 = np = 10 ^ {15}),勉强可以接受。

模板题 P4245 任意模数多项式乘法代码。注意用 long doubledouble 会被卡精度,或者换一种精度较高的 FFT 写法

5.2.2 三模数 NTT

前置知识:(扩展)中国剩余定理。

选取三个常用 NTT 模数,分别算出 (f * g) 在这些模数下的结果,再使用中国剩余定理合并即可。

我们选择 (p_1 = 998244353)(p_2 = 1004535809)(p_3 = 469762049),设结果分别为 (h_1, h_2, h_3)

如果使用 CRT,则需要 __int128,因为 (h) 的真实值为

[(h_1p_2p_3times mathrm{inv}(p_2p_3, p_1) + h_2p_1p_3times mathrm{inv}(p_1p_3, p_2) + h_3p_1p_2times mathrm{inv}(p_1p_2, p_3)) bmod (p_1p_2p_3)
]

使用 exCET 就不需要 __int128

  • 合并 (hequiv h_1pmod {p_1})(hequiv h_2pmod {p_2})。设 (h = h_1 + yp_1),则 (h_1 + yp_1 equiv h_2pmod {p_2}),解得 (yin [0, p_2)) 之后得到 (h equiv h_1 + yp_1 pmod {p_1p_2}),记作 (hequiv xpmod {p_1p_2})
  • 合并 (h equiv xpmod {p_1p_2})(hequiv h_3 pmod {p_3})。设 (h = x + yp_1p_2),类似解得 (yin [0, p_3)),得到 (h equiv x + yp_1p_2pmod {p_1p_2p_3})。因为 (x + yp_1p_2 ,所以它就是真实值,可以直接取模。

效率比拆系数 FFT 低不少,因为进行了九次 DFT。代码

5.3 分治 NTT

前置知识:CDQ 分治。

分治 NTT 在信息竞赛界应用广泛,这归功于他所解决的问题:形如 (f_i = sum_{j = 0} ^ {i - 1} f_j g_{i - j})半在线卷积

5.3.1 算法介绍

形如给定 (g_1sim g_{n - 1})(f_0),求 (f_1sim f_{n - 1}) 满足 (f_i = sum_{j = 0} ^ {i - 1} f_j g_{i - j}) 的问题被称为 半在线卷积。因为 (f_i) 很大,一般会对 NTT 模数 (p) 取模。

我们不能通过简单的 NTT 解决这个问题,因为 (f) 的每一项均和之前项有关,这要求我们在尝试计算 (f_i) 时必须已经求出 (f_0sim f_{i - 1}),而 NTT 做不到这一点。或者说,它们解决的问题形式不同。

这是一个在线的问题,考虑使用 CDQ 分治 转化为静态问题。

  • 设当前分治区间为 ([l, r]),其中 (f_0 sim f_{l - 1})(f_lsim f_r) 的贡献已经计算完毕。
  • (l = r),说明已经求出 (f_l),返回。
  • 否则,令 (m = lfloor frac {l + r} 2rfloor),先向左递归 ([l, m]) 求解 (f_lsim f_m)
  • 为了求解 (f_{m + 1}sim f_r),我们需要计算 (f_lsim f_m) 对它们的贡献:(f_i = sum_{j = l} ^ m f_j g_{i - j})。因为 (f_lsim f_m) 已知,所以总的贡献相当于两个已知的多项式相乘的结果。具体地,令 (f'_j = f_{j + l}(0leq j leq m - l)),计算 (h = f' * g),则 (f_i) 受到 (h_{i - l}) 的贡献。
  • 向右递归 ([m + 1, r]) 求解 (f_{m + 1}sim f_r)

因为 (f, g) 的长度均不超过区间长度,所以时间复杂度 (mathcal{O}(nlog ^ 2 n))

模板题 P4721 分治 FFT 代码:

#include 
using namespace std;

using ll = long long;
using ull = unsigned long long;

constexpr int N = 1 = mod && (x -= mod);}
int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % mod;
    a = 1ll * a * a % mod, b >>= 1;
  }
  return s;
}

int n, f[N], g[N];
void solve(int l, int r) {
  if(l == r) return;
  int m = l + r >> 1, L = 1;
  solve(l, m);
  while(L > n;
  for(int i = 1; i 

5.3.2 阶梯网格路径计数

阶梯网格路径计数是一类可以使用分治 NTT 解决的经典问题。

给定 (n + 1)(m + 1) 行的网格图,左下角和右上角分别记为 ((0, 0))((n, m))。从 ((0, 0)) 出发,每次只能向右或向上走,求走到 ((n, m)) 的方案数。让我们回忆:方案数为 (binom {n + m} n)

当然,问题没有这么简单。我们限制第 (i) 列不能走到行数大于 (c_i) 的点,其中 (c_i) 单调不降,且 (c_n = m)。换言之,求出在一个阶梯网格从左下角走到右上角的方案数。

显然有 (mathcal{O}(nm)) 的动态规划:设 (f_{i, j}) 表示走到 ((i, j)) 的方案数,则 (f_{i, j}) 转移至 (f_{i + 1, j})(f_{i, j + 1})

考虑一段 (c_i) 相同的极长区间 ([l, r](l 从 (f_{l, j_1}) 转移到 (f_{r, j_2}) 的贡献系数:为防止重复计数,从 ((l, j_1)) 出发的第一步应当向右,因此有系数 (binom {(r - l - 1) + (j_2 - j_1)} {r - l - 1})。设 (g_i) 表示 (binom {r - l - 1 + i} {i}),则 (f_{l, j_1} g_ito f_{r, j_1 + i}),为卷积形式。

在不受 (c_i) 影响的时候,我们确实可以这样做。但是对于每个 (i),它内层的所有 (j) 之间也有转移,这让我们不方便借助分治和卷积加快整个过程。考虑进行一些等价转换便于分层。

稍作思考,设 (f_{k, j}) 表示走到 ((i, j)) 其中 (i + j = k) 的方案数,则 (f_{k, j}) 转移至 (f_{k + 1, j})(f_{k + 1, j + 1})

进一步地,为了避免在 (k > n) 时受到 (jgeq k - n) 的限制,考虑从 ((n, 0)) 开始沿左上方向把整个阶梯网格劈成两半,分别计算从 ((0, 0))((n, m)) 到斜对角线上的点((i + j = n))的方案数,而这两个问题是对称的。

综合上述分析,我们将问题转化为:从 ((0, 0)) 出发,每次可以向右或右上走一步,求走到最右侧一列每个点的方案数。其中,在第 (i) 列不能走到行数大于 (c_i) 的点。这个 (c_i) 可通过原问题的 (c_i) 进行简单转化得到,且容易证明其仍不降且满足很好的性质:(c_{i + 1} - c_i leq 1)

因此,从 (f_l) 推到 (f_r) 的时候,我们会发现对于 (jleq c_r - (r - l)),设不等号右侧的数为 (x)(f_{l, j})(f_r) 的贡献不会受到 (c) 的影响:因为 (c_{i + 1} - c_ileq 1),所以从 ((l, j)) 开始,就算每一步都往右上走,也不会突破限制。这样,(f_{l, 0}sim f_{l, x})(f_r) 的贡献可直接卷积计算。对于 (f_{l, x + 1}sim f_{l, c_l}),分治下放至左右子区间 ((l, m])((m, r]) 进行递归。

多项式 I:拉格朗日插值与快速傅里叶变换插图2

有了大致思路,设计算法就很简单了:设 (mathrm{solve}(l, r, Delta, F)) 表示当前区间为 ((l, r]),传入多项式 (F) 的第 (i) 项表示 (f_{l, i + Delta}),返回多项式 (G) 的第 (i) 项表示 (f_{r, i + Delta})。也可视作暂时将 (c_lsim c_r) 全部减去 (Delta),传入 (f_l = F),返回 (G = f_r),接下来就使用这种视角。

对于 (j leq c_r - (r - l)),设不等号右侧的数为 (x)(F_0sim F_x)(H_0sim H_{r - l}) 卷积求出 (F_0sim F_x) 转移得到的 (G_0sim G_{c_r}),其中 (H_i) 即转移系数 (binom {r - l} i)

对于 (j > x),分治下放:设 (F' = F_{x + 1}sim F_{c_l})(mid = lfloor frac {l + r} {2} rfloor),则先递归左半部分 (F'gets mathrm{solve}(l, mid, Delta + x + 1, F')),再递归右半部分 (F'gets mathrm{solve}(mid, r, Delta + x + 1, F')),则此时 (F') 表示从 (F_{x + 1}sim F_{c_l}) 转移得到的 (G_{x + 1}sim G_{c_r})

两部分相加即得 (G),返回即可。

边界 (r - l = 1) 的处理是平凡的。

两次下传 (F') 时,(F') 的长度显然不超过 (c_r - x = r - l),因此在处理区间 ((l, r]) 时涉及到的所有多项式的长度均不超过其父区间长度。设 (n, m) 同级,则时间复杂度为 (mathcal{O}(nlog ^ 2 n))

接下来对问题进行一些扩展:

  • 如果没有 (c_{i + 1} - c_ileq 1) 的性质,但 (c_i) 单调不降,还能做吗?答案是可以,只需将 (x) 的定义改为 (c_l - (r - l)),此时 ((l, r]) 涉及到的多项式长度不超过父区间的长度加上端点处 (c) 的差值,可以接受。

5.4 例题

CF1770G Koxia and Bracket

相比于 F,这道题就略显套路了。

将左括号视为 (1),右括号视为 (-1),找到最后一个使得前缀和小于 (0) 的位置 (lst),则 (s_1sim s_{lst}) 的每个左括号都有用,必然不会删去。同理,(s_{lst + 1}sim s_n) 的每个右括号也都有用。

对于 (s_1sim s_{lst}) 的每个右括号 (s_i),如果它对应的前缀和为 (c_i),则为保证前缀和非负,在 (s_1sim s_i) 中至少需要删掉 (-c_i) 个右括号。我们只关心 (-c_i) 的前缀最大值,因为若这些位置满足条件,则所有位置满足条件,而每个前缀最大值一定会比先前的前缀最大值大 (1),所以我们将问题转化为:给定长为 (m) 的序列,其中有 (c) 个位置被打上了标记,求出选择 (c) 个位置的方案数,使得对于每个前缀,选择的位置数不小于被打上标记的位置数。

(f_{i, j}) 表示考虑到第 (i) 个位置,当前选择位置数减去打标记位置数的数量为 (j)。对于没有被打标记的位置 (p)(f_{p - 1, j}) 转移到 (f_{p, j / j + 1}),否则转移到 (f_{p, j - 1 / j})

非常明显的阶梯格路计数问题,稍有变形,不过解法大差不差,核心思想是一致的:考虑 (f_l) 转移到 (f_r),用卷积描述一部分转移,剩下来 (mathcal{O}(r - l)) 个位置分别递归两侧处理。

对于本题,就是 (f_{l, j}(jgeq cnt))(f_r) 的贡献用卷积描述,其中 (cnt) 表示 (l + 1sim r) 的被打上标记的位置。剩余部分递归,长度只有 (cnt),而且因为截取的是低位,我们甚至不需要记录偏移量 (Delta)

(r - l = 1) 的边界情况是平凡的。

(s_{lst + 1}sim s_n) 的左括号类似处理即可。

每个分治区间涉及到的多项式长度不超过其父区间长度,因此时间复杂度为 (mathcal{O}(nlog ^ 2n))代码

参考资料

第一章:

第二章:

第四章:

第五章:

        
posted @

2023-01-08 18:23 
qAlex_Weiq 
阅读(274
评论(5
编辑 
收藏 
举报

文章来源于互联网:多项式 I:拉格朗日插值与快速傅里叶变换

THE END
分享
二维码