[多项式学习笔记] 拉格朗日插值

[多项式学习笔记] 拉格朗日插值

多项式插值

给定 \(x\) 坐标两两不同的 \(n + 1\) 个点，能够唯一确定一个 \(n\) 次多项式。从给定点求出多项式的过程称为插值。

具体而言，给定 \(n + 1\) 个点 \((x_0, y_0), (x_1, y_1), \cdots, (x_n, y_n)\)，求 \(n\) 次多项式满足

\[f(x_i) = y_i, \forall 0 \le i \le n \]

拉格朗日插值法可以在 \(O(n^2)\) 的时间内求出这个多项式。

考虑构造 \(n + 1\) 个 \(n\) 次多项式 \(f_0(x), f_1(x), \cdots, f_n(x)\)，第 \(i\) 个多项式满足

\[f_i(x_j) = \begin{cases} y_j, \ j = i \\ 0, \ \text{Otherwise} \end{cases} \]

也就是它的图像经过第 \(i\) 个点，和其他点在 \(x\) 轴上的投影。易知待求函数 \(f(x) = \sum_{i = 1}^{n} f_i(x)\)。

由于 \(f_i(x)\) 经过除了 \(x_i\) 的其它点在 \(x\) 轴上的投影，所以可以设 \(f_i(x) = a \cdot \prod_{j \neq i} (x - x_j)\)。根据 \(f_i(x_i) = y_i\)，解得 \(a = \dfrac{y_i}{\prod_{j \neq i} (x_i - x_j)}\)。因此

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

从而

\[f(x) = \sum_{i = 0}^{n} y_i \prod_{j \neq i}\dfrac{x - x_j}{x_i - x_j} \]

于是就可以在 \(O(n^{2})\) 的时间内暴力求出 \(f(x)\) 当 \(x = k\) 时的点值，其中 \(k\) 是任意的数。注意，我们现在给出了 \(f(x)\) 的一个表示，但这并不是 \(f(x)\) 的系数表示法。而如果能求出 \(f(x)\) 的系数表示法，就可以在 \(O(n)\) 的时间内求点值。在需要多次计算点值时，这是很有用的。下面将介绍如何求出它。

插值求出系数表示法

由于 \(f(x)\) 是 \(n\) 次多项式，所以可以表示为

\[f(x) = \sum_{i = 0}^{n} a_{i} x^{i} \]

我们希望求出 \(a_{0..n}\)。

显然，只需对每个 \(f_i(x)\) 求出其展开后 \(x^{k}\) 的系数就可以了。把 \(f_{i}(x)\) 改写成

\[f_i(x) = \dfrac{y_{i}}{\prod_{j \neq i} (x_i - x_j)} \times \prod_{j \neq i} (x - x_j) \]

前一项是常数，容易求出。我们主要关心第二项。

记 \(F(x) = \prod_{j = 1}^{n} (x - x_j)\)。我们先求出 \(F(x)\) 的系数再对每个 \(i\) 求出 \(f_i(x)\) 中第二项的系数。

直接展开是困难的，我们使用递推解决这个问题。记 \(f(i, k)\) 表示 \(\prod_{j = 1}^{i} (x - x_{j})\) 展开以后 \(x^{k}\) 的系数。初始化 \(f(0, 0) = 1\)，转移方程为

\[f(i, k) = f(i - 1, k - 1) - x_{i} \cdot f(i - 1, k) \]

这是因为 \(x^{k}\) 的系数有两个来源：一是 \(x\) 乘上 \(\prod_{j = 1}^{i - 1} (x - x_{j})\) 中的 \(k - 1\) 次项，这不改变系数；而是 \(-x_i\) 乘上 \(\prod_{j = 1}^{i - 1} (x - x_{j})\) 中的 \(x^{k}\) 项，这把系数变为原来的 \(-x_i\) 倍。

于是我们就可以在 \(O(n^2)\) 的时间内求出 \(F(x)\) 展开以后各项的系数。现在我们还要对每个 \(i\) 求出 \(\dfrac{F(x)}{x - x_i}\) 展开以后的系数。根据定义，\(\dfrac{F(x)}{x - x_i}\) 一定也是整式。记 \(F(x) = (x - x_i)P(x)\)，我们就是要求出 \(P(x)\) 的系数表示。记 \(b_k = [x^{k}]P(x)\)，把多项式的竖式乘法写出来：

我们得到了 \(-b_0x_i = a_0\)，于是就可以解出 \(b_0\)。对于 \(1 \le j \le n\) 又有 \(a_{j - 1} = b_{j - 1} - b_{j}x_{i}\)，因此我们可以递推求出 \(b_{j}\)，这样就求出了 \(P(x)\) 的系数表示，进而求出了 \(f_i(x)\) 的和 \(f(x)\) 系数表示，总时间复杂度为 \(O(n^{2})\)。

横坐标是连续整数的插值

如果横坐标是连续整数，我们可以做到 \(O(n)\) 插值。

设已知 \(f(1), f(2), \cdots, f(n + 1)\)，代入插值公式得

\[f_{i} = y_{i} \prod_{j \neq i} \dfrac{x - j}{i - j} \]

分别处理分子和分母。分子为

\[\prod_{j \neq i} (x - j) = \dfrac{\prod_{j} (x - j)}{x - i} \]

（注意，\(x - i\) 出现在了分母中，这意味着如果我们想求的点值和某个 \(i\) 相同时需要特判）

分母为

\[\begin{aligned} \prod_{j \neq i} (i - j) &= \prod_{1 \le j < i} (i - j) \times \prod_{i < j \le n + 1} (i - j) \\ &= (-1)^{n + 1 - i}(i - 1)!(n + 1 - i)! \end{aligned} \]

因此

\[f_{i}(x) = (-1)^{n + 1 - i} y_{i} \dfrac{\prod_{1 \le j \le n + 1} (x - j)}{(x - i)(i - 1)!(n + 1 - i)!} \]

预处理阶乘和阶乘逆元后就可以在 \(O(1)\) 的时间内计算每个 \(f_i(x)\)，从而在 \(O(n)\) 的时间内计算 \(f(x)\)。

例题

I. CF622F. The Sum of the k-th Powers

（题目链接）求 \(\sum_{i = 1}^{n} i^{k}\)，对 \(10^{9} + 7\) 取模。

\(1 \le n \le 10^{9}\)，\(0 \le k \le 10^{6}\)。

本题是经典的自然数幂和问题，解决的方法有很多。下面使用的是《具体数学》中介绍的扰动法。

记 \(S_{k}(n) = \sum_{i = 0}^{n} i^{k}\)，即自然数的 \(k\) 次方和。那么

\[\begin{aligned} S_{k + 1}(n) + (n + 1)^{k + 1} &= \sum_{0 \le i \le n + 1} i^{k + 1} \\ &= \sum_{1 \le i \le n} (i + 1)^{k + 1} \\ &= \sum_{0 \le i \le n} \sum_{0 \le j \le k + 1} \dbinom{k + 1}{j}i^{j} \quad \text{（二项式定理）} \\ &= \sum_{0 \le j \le k + 1} \dbinom{k + 1}{j} \sum_{0 \le i \le n} i^{j} \quad \text{（交换求和顺序）} \\ &= \sum_{0 \le j \le k + 1} \dbinom{k + 1}{j} S_{j}(n) \\ &= S_{k + 1}(n) + (k + 1)S_{k}(n) + \sum_{0 \le j < k} \binom{k + 1}{j} S_{j}(n) \end{aligned} \]

因此

\[\boxed{S_{k}(n) = \dfrac{(n + 1)^{k + 1} - {\large\sum_{j = 0}^{k - 1}} \binom{k + 1}{j} S_{j}(n)}{k + 1}} \]

验证一下，当 \(k = 2\) 时有

\[\begin{aligned} \sum_{i = 0}^{n} i^{2} &= \dfrac{(n + 1)^{3} - (S_{0}(n) + 3S_{1}(n))}{3} \\ &= \dfrac{n^{3} + 3n^{2} + 3n + 1 - (n + 1 + \dfrac{3(n + 1)n}{2})}{3} \\ &= \dfrac{n(n + 1)(2n + 1)}{6} \end{aligned} \]

这给出了我们熟悉的结果。

由于 \((n + 1)^{k + 1}\) 的存在，所以 \(S_{k}(n)\) 至少是 \(n\) 的 \(k + 1\) 次多项式。如果 \(S_{j}\)（\(0 \le j \le k - 1\)）没有贡献 \(n\) 的更高次项，就可以确定：\(S_k(n)\) 是 \(n\) 的 \(k + 1\) 次多项式。

可以用数学归纳法证明这一点：因为 \(S_{0}(n) = n + 1\) 是关于 \(n\) 的 \(1\) 次多项式，所以基础成立，进而上述结论成立。

得到这个结论之后，我们就可以用拉格朗日插值求出这个多项式，并计算它的点值。

不妨选择计算 \(1 \sim k + 2\) 的点值，这样就可以套用连续整数的拉插。由于 \(f(x) = x^{k}\) 是（完全）积性函数，可以用线性筛在 \(O(k)\) 的时间内计算点值。因此总时间复杂度为 \(O(k \log P)\)，其中 \(P\) 是模数。AC 记录

II. [COTS 2023] Mapa

题目链接：洛谷 | QOJ

我们只有 \(3000\) 个 bit 可用，这最多只能表示 \(100\) 个 \(10^{9}\) 以内的非负整数，似乎不可能记录下 \(100\) 个键值对。解决方法非常巧妙：我们构造 \(n - 1\) 次多项式函数 \(f(x)\)，满足对于每个 \(1 \le i \le n\) 都有 \(f(x_i) = y_i\)。我们只需要 \(n\) 个点值就可以还原这个多项式，进而可以求出每个键对应的值。不妨记录 \(x = 1, 2, \dots, n\) 时的点值，这样就可以存下了。

因此，编码时，对于 \(1 \le x \le n\)，使用拉插在 \(O(n^{2})\) （求逆元的时间不计，下同）的时间内求出 \(f(x)\) 的值，总时间复杂度 \(O(n^3)\)（。解码时，对于每个询问的 key \(x\)，使用拉插在 \(O(n^{2})\) 的时间内求出 \(f(x)\)，这也就对应了要求的 value。

AC 记录

一些数学细节：

1. 为了避免溢出，计算的过程都在模一个大于 \(10^{9}\) 的质数的意义下进行。只要模数大于 \(10^{9}\)，我们就不会损失信息。
2. 可以证明当 \(x\) 是整数时 \(f(x)\) 都是整数。（？）

III. [CrCPC 2024] EA Enigma

（题目链接）

由于我们一开始对序列完全不了解，所以可以随便猜一个序列，不妨猜一个全为 \(1\) 的序列。假设我们猜中了 \(i\) 个位置，那么还剩下 \(n - i\) 个位置不确定，但我们成功把这些位置可能的值的个数减少了 \(1\)，于是递归到序列长度为 \(n - i\)，值域为 \(k - 1\) 的子问题。设 \(F(n, k)\) 表示序列长度为 \(n\)，值域为 \(k\) 时的答案，则有如下递归式：

\[F(n, k) = 1 + \sum_{i = 0}^{n} \dfrac{\binom{n}{i}(k - 1)^{n - i}}{k^{n}} F(n - i, k - 1) \]

计算这个式子的时间复杂度为 \(O(nk)\)，不过它确实是正确的，因为我写了发 暴力 验证了一下。但它的形式过于复杂，还是递归式，难以化简，所以根本就没有前途（悲）。

在一开始的时候就不要列这种式子。沿着一开始的思路：如果第一次猜测以后，还有位置没有确定，不妨猜测剩下的位置都为 \(2\)。实际上只要不猜重复，猜什么都是一样的。然后再猜 \(3\)，\(4\)，一直猜到完全确定，于是猜测的次数就是序列的最大值。也就是说猜一个序列的次数等于序列中的最大值。

根据经典的容斥，长度为 \(n\) 且最大值为 \(i\) 的序列有 \(i^{n} - (i - 1)^{n}\) 个，因此

\[ans = \dfrac{\sum_{i = 1}^{k} i(i^{n} - (i - 1)^{n})}{k^{n}} \]

式子中有对 \(k\) 的求和，必须去掉：

\[\begin{aligned} k^{n}ans &= \sum_{i = 1}^{k} i(i^{n} - (i - 1)^{n}) \\ &= \sum_{i = 1}^{k} i \left(i^{n} - \sum_{j = 0}^{n} \dbinom{n}{j} (-1)^{n - j}i^{j} \right) \\ &= \sum_{i = 1}^{k} i \left(\sum_{j = 0}^{n - 1} \dbinom{n}{j} (-1)^{n - j + 1}i^{j} \right) \\ &= \sum_{j = 0}^{n - 1} (-1)^{n - j + 1} \dbinom{n}{j}\sum_{i = 1}^{k} i^{j + 1} \\ &= \sum_{j = 0}^{n - 1} (-1)^{n - j + 1} \dbinom{n}{j} S_{j + 1}(k) \\ &= \sum_{j = 0}^{n - 1} (-1)^{n - j + 1} \dbinom{n}{j} O(k^{j + 2}) \\ &= O(k^{n + 1}) \end{aligned} \]

推这个式子的过程中，一开始看到了对 \(i(i^{n} - (i - 1)^{n})\) 的求和，这看上去和自然数幂和的形式很相似，于是我们猜测 \(k^{n}ans\) 是一个关于 \(k\) 的大约 \(n + 1\) 次的多项式，最终的结果确实如此。然后就可以用拉格朗日插值求出答案。总时间复杂度为 \(O(n (\log n + \log P))\)，其中 \(P\) 是模数。（如果用线性筛，就可以去掉 \(\log n\)。）AC 记录

IV. [集训队互测 2012] calc

题目链接

先考虑怎么设计状态。为了方便，我们只研究递增序列（或者说把序列看作元素的集合，本质上都是把有序转化成无序）的方案数，计算答案时再乘上 \(n!\)。

设 \(f(i, j)\) 表示长度为 \(i\) 且最大值不超过 \(j\) 的合法序列数，分最大值为 \(j\) 和最大值小于 \(j\) 两种情况得到转移方程：

\[f(i, j) = j \cdot f(i - 1, j -1 ) + f(i, j -1 ) \]

初始化 \(f(0, 0) = 1\)。

则答案为 \(f(n, k) \cdot n!\)，直接计算的时间复杂度为 \(O(nk)\)。