拉格朗日插值
本文作者为 JustinRochester。
拉格朗日插值
拉格朗日插值的基础
拉格朗日插值旨在给定 \(n+1\) 个互不相同的变量值 \(x_i(i=0,1,2,\cdots, n)\) 与他们对应的 \(y_i\) 的情况下，构造出一个至 \(n\) 次的多项式函数 \(f(x)\) ，满足 \(f(x_i)=y_i\) 。
这个算法的妙用是能解决一类多项式函数的数列题。
如给定数列 \(a_1=1, a_2=2, a_3=4, a_4=8\) 问数列第 \(5\) 项是什么？很显然是 \(a_5=114514\) ，因为他们分别是多项式函数 \(\displaystyle f(x)={114499\over 24}x^4-{190831\over 4}x^3+{4007453\over 24}x^2-{954153\over 4}x+114499\) 在 \(x\) 分别取 \(1,2,3,4,5\) 时的函数值。
额，等一下，这个并不是妙用。
如 Codeforces Round #841 (Div. 2) and Divide by Zero 2022 的 B. Kill Demodogs ，它打表找到规律后发现是 \(\displaystyle \sum_{i=1}^n (i(i-1)+i^2)\) 。
如果你并不会自然数平方和的公式或者有限微积分科技，你可能并不能求解出这个答案。
但根据分析，二次方求和的公式是至多三次的多项式可以描述的，因此你可以暴力打出前 \(4\) 项结果，然后采用拉格朗日插值的方法通过本题。
构造“全 \(0\) ”多项式
我们不妨先抛开整个拉格朗日插值的要求，先考虑一个简单的问题：
给定 \(n\) 个互不相同的变量值 \(x_i\) ，如何构造一个多项式函数 \(f(x)\) ，使得 \(f(x_i)\) 均为 \(0\) 呢?
根据中学知识，我们很显然可以构造出多项式 \(\displaystyle f(x)=\prod_i (x-x_i)\) 。这显然满足要求。
当然，任意满足 \(\displaystyle f(x)=\omega(x)\prod_i (x-x_i)\) 的多项式函数均满足要求，其中 \(\omega(x)\) 是任意多项式函数。
此外，我们可以稍微分析一下，得知 \(f(x)\) 是一个至多 \(n\) 次多项式。
根据代数基本定理，任意 \(n\) 次多项式至多只有 \(n\) 个复数根。因此，\(f(x)\) 仅有这 \(n\) 个自变量取值处，函数值为 \(0\) ，否则函数值均不为 \(0\) 。
构造基多项式
我们现在在“全 \(0\) ”多项式的要求上进一步拓展，现在我们要求：
给定 \(n+1\) 个互不相同的变量值 \(x_i\) ，如何构造一个多项式函数 \(l_p(x)\) ，使得 \(l_p(x_i)\) 在 \(i\neq p\) 时均为 \(0\) ，而 \(i=p\) 时为 \(1\) 呢？
根据上文，我们可以很快的构造出一个多项式函数满足第一个要求： \(\displaystyle f_p(x)=\prod_{i\neq p}(x-x_i)\) 。
但很可惜，它只能保证 \(f_p(x_p)\neq 0\) ，并不能保证 \(f_p(x_p)=1\) 。
但是，我们可以意识到，只需要将 \(f_p(x)\) 函数整体除以 \(f_p(x_p)\) 即可得到答案。因为 \(f_p(x_p)\neq 0\) ，因此它是可以成为分母的。
于是我们构造出了基多项式 \(\displaystyle l_p(x)={f_p(x)\over f_p(x_p)}\) 。
这里，我们也可以将 \(f_p(x_p)\) 按 \(f_p(x)\) 的定义展开，我们就能得到： \(\displaystyle l_p(x)=\prod_{i\neq p}{x-x_i\over x_p-x_i}\)
此外，我们可以在这基础上进行任意的拓展：
给定 \(n+1\) 个互不相同的变量值 \(x_i\) 和常数 \(y_p\)，如何构造一个多项式函数 \(l_p(x)\) ，使得 \(f(x_i)\) 在 \(i\neq p\) 时均为 \(0\) ，而 \(i=p\) 时为 \(y_p\) 呢？
那我们可以直接写下一个答案： \(l_p(x)\cdot y_p\) 。
拉格朗日插值的公式
好，经过上述两步我们重新回到原问题：
给定 \(n+1\) 个互不相同的变量值 \(x_i(i=0,1,2,\cdots, n)\) 与他们对应的 \(y_i\) 的情况下，如何构造出一个至 \(n\) 次的多项式函数 \(f(x)\) ，满足 \(f(x_i)=y_i\) 呢？
这个问题并不好直接解决，但如果我能分别构造只在 \(x_0\) 处取值为 \(y_0\) （其他位置取值为 \(0\) ，下同）、只在 \(x_1\) 处取值为 \(y_1\) 、...、只在 \(x_n\) 处取值为 \(y_n\) 的多项式。我们分别求和，就能得到答案了。
于是，根据上面基多项式的拓展，我们可以直接写下了拉格朗日插值得到答案： \(\displaystyle L(x)=\sum_i y_il_i(x)\) 。
显然，我代入任何一个 \(x_p\) 时，\(\displaystyle L(x_p)=\sum_i y_il_i(x_p)=0+y_pl_p(x)+0=y_p\cdot 1=y_p\) 。
同理，我们也可以将它展开，即可得到 \(\displaystyle L(x)=\sum_i y_i\prod_{j\neq i}{x-x_j\over x_i-x_j}\) 。
运用这个方法，我们就可以根据 \(n+1\) 个点，拟合出满足它们的多项式函数了。
此外，在自变量取值点连续，为 \(t, t+1, \cdots t+n\) 时，拉格朗日插值函数可进一步进行简化：
\(\begin{aligned} &L(x) \\=&\sum_i y_i\prod_{j\neq i}{x-x_j\over x_i-x_j} \\=&\sum_i y_i\prod_{j\neq i}{x-(t+j)\over (t+i)-(t+j)} \\=&\sum_i y_i\prod_{j\neq i}{x-t-j\over i-j} \\=&\sum_i y_i\prod_{j=0}^{i-1}{1\over i-j}\prod_{j=i+1}^n{-1\over j-i}\prod_{j\neq i}(x-t-j) \\=&\sum_i {y_i\cdot (-1)^{n-i}\over i!(n-i)!}\prod_{j\neq i}(x-t-j) \end{aligned}\)
特别地，当 \(x\not\in\{t, t+1, \cdots, t+n\}\) 时，\(\displaystyle \prod_{j\neq i}(x-t-j)={1\over x-t-i}\prod_j (x-t-j)={(x-t)^{\underline {n+1}}\over x-t-i}\)
其中，\(x^{\underline n}\) 表示 \(x\) 的 \(n\) 次下降幂，即 \(x^{\underline n}=x(x-1)(x-2)\cdots (x-n+1)\) 。
于是此时，\(\displaystyle L(x)=\sum_i {y_i\cdot (-1)^{n-i} (x-t)^{\underline {n+1}}\over i!(n-i)!(x-t-i)}\) 。
尤其是当我们自己打表时，我们会尽可能打出 \(0\) 到 \(n\) 的取值，然后插值后面的答案。因此，该公式退化为 \(t=0\) 的情况 \(\displaystyle L(x)=\sum_i {y_i\cdot (-1)^{n-i} x^{\underline {n+1}}\over i!(n-i)!(x-i)}\) 。
拉格朗日插值的唯一性
这一节并不是说插值出来的函数是唯一的。实际上，根据“全 \(0\) ”多项式，我们仍可以插出在这 \(n+1\) 个点上均为 \(0\) 的多项式（并且有无穷多个）。我们挑选任何一个和上一节讲到的拉格朗日插值函数相加，都是满足题设的。
这一节说的是，在不超过 \(n\) 次的多项式中，这个插值的结果是唯一的。
我们不妨假设这不超过 \(n\) 次的多项式是 \(\displaystyle P(x)=\sum_{i=0}^n a_ix^i\) 。其中任意的 \(a_i\) 均可能为 \(0\) 。
我们在分别求值 \(P(x_i)=y_i\) 的过程，实际上可以用矩阵的形式进行表现：
\( \begin{aligned} \left( \begin{matrix} y_0\\ y_1\\ \vdots\\ y_n\\ \end{matrix} \right) =\left( \begin{matrix} x_0^0 & x_0^1 & \cdots & x_0^n\\ x_1^0 & x_1^1 & \cdots & x_1^n\\ \vdots & \vdots & \ddots & \vdots\\ x_n^0 & x_n^1 & \cdots & x_n^n\\ \end{matrix} \right) \cdot \left( \begin{matrix} a_0\\ a_1\\ \vdots\\ a_n \end{matrix} \right) \end{aligned} \)
我们可以把矩阵记为 \(Y=XA\) 的形式。
而注意到矩阵 \(X\) ，它是一个标准的 Vandermonde 矩阵，行列式 \(\displaystyle \det (X)=\prod_{1\leq i<j\leq n} (x_j-x_i)\) 。
由于我们假设了自变量的取值互不相同，故 \(x_j-x_i\neq 0\) ，因此 \(\det (X)\neq 0\) ，矩阵 \(X\) 可逆。
于是我们已知矩阵 \(Y, X\) ，求矩阵 \(A\) ，则 \(A=X^{-1}Y\) ，这个 \(A\) 矩阵必然是唯一的。
因此可以肯定，拉格朗日插值得到的函数在不超过 \(n\) 次的多项式函数中，是唯一的。
这也保证了我们分析后，插值出来的函数必定是我们要求的函数。
如何优雅的插值
以下介绍一些拉格朗日插值的做题技巧。
首先考虑到原拉格朗日插值函数 \(\displaystyle L(x)=\sum_i y_i\prod_{j\neq i}{x-x_j\over x_i-x_j}=\sum_i {y_i\over \prod_{j\neq i}(x_i-x_j)}\prod_{j\neq i}(x-x_j)\) 。
由于可以通过构造，使得 \(x_i-x_j\) 的取值有 \(O(n^2)\) 种。因此该算法无论如何都无法在 \(O(n)\) 的时间内插值出某个点的结果。
但是，我们可以在 \(O(n^2)\) 时间内计算出 \(\displaystyle m_i=\prod_{j\neq i}(x_i-x_j)\) 。然后我们想要在尽可能快的时间内求序列 \(m_i\) 每一项的逆元，这可以用线形求逆元的办法：在 \(O(n)\) 的时间内计算出 \(m_i\) 的前缀积，花 \(O(\log n)\) 的时间计算出 \(m_i\) 乘积的倒数，最后再 \(O(n)\) 向前递推，算出 \({1\over m_i}\) ，并计算出 \({y_i\over m_i}\) 。
预处理的时间复杂度为 \(O(n^2)+(O(n)+O(\log n)+O(n))=O(n^2)+O(n)=O(n^2)\) 。
而后续查询 \(L(x_k)\) 时，我们可以先查询是否 \(x_k\) 为插值点，是则直接输出答案。否则我们可以通过 \(O(n)\) 预处理 \(x_k-x_j\) 数组的前缀积 \(pre_k\) 和后缀积 \(suf_k\)，然后直接 \(O(n)\) 计算答案 \(\displaystyle L(x_k)=\sum_i {y_i\over m_i}\cdot pre_{i-1}suf_{i+1}\) 。
当查询 \(T\) 次时，总复杂度为 \(O(n^2+Tn)\) 。
而对于自己打表或连续插值时，我们会采用公式 \(\displaystyle L(x)=\sum_i {y_i\cdot (-1)^{n-i} (x-t)^{\underline {n+1}}\over i!(n-i)!(x-t-i)}\) 。
同理，当查询的值为插值点时，直接查表输出答案，否则对于 \({y_i\cdot (-1)^{n-i}\over i!(n-i)!}\) 可以用上述办法在 \(O(n)\) 内预处理；\((x-t)^{\underline {n+1}}\) 也可以在 \(O(n)\) 时间内计算，而 \(1\over x-t-i\) 即要求我们线性求逆元，也可以用上文的方法在 \(O(n)\) 内求解。
因此，当查询 \(T\) 次时，总复杂度为 \(O(n+Tn)=O(Tn)\) 。
