【学习笔记】拉格朗日插值
基本内容
考虑一个最高次项为\(n\)的多项式函数\(f(x)\)，它可以写成：
拉格朗日插值要解决的问题是：
告诉你\(n+1\)个\(x_i\)(\(0\leq i\leq n\))，和它们对应的\(y_i=f(x_i)\)的值。让你回答询问：对某个给定的\(k\)，求出\(f(k)\)的值。
朴素的想法是把\(a_0,a_1\dots a_n\)看做\(n+1\)个未知数，给定的\(n+1\)组\(x_i\)和\(y_i\)看做\(n+1\)个方程。然后用高斯消元解方程。时间复杂度是\(O(n^3)\)。
考虑更快的做法。也就是拉格朗日插值法。它的结论是：
可以感性理解其正确性：当\(k\)等于某个\(x_t\)时：若\(i\neq t\)，则必会出现\(j=t\)的情况，也就是说\(\prod\)中必有一项的分子为\(0\)，因此\(i\neq t\)的项对求和式就没有贡献。当\(i=t\)时，\(\prod\)的每一项分子分母都相等。因此\(f(k)=y_t\)，这与题目条件是相符的。
至此，我们已经求出了\(f(k)\)。时间复杂度就是按这个式子暴力模拟的复杂度：\(O(n^2)\)。
求系数
我们不满足于求某个特定\(f(k)\)的值。我们还想要知道\(a_0\dots a_n\)具体是多少。
把上面的式子转化一下：\(f(k)=\sum_{i=0}^{n}\frac{y_i}{\prod_{j\neq i}x_i-x_j}\prod_{j\neq i}(k-x_j)\)。请务必注意：这个式子里只有\(k\)是变量，其他值都是常量。也就是说，我们的目标是要把它写成\(f(k)=a_0+a_1k+a_2k^2+\dots a_nk^n\)的这种形式。
其中\(\frac{y_i}{\prod_{j\neq i}x_i-x_j}\)，作为常量，可以直接求出来。
那么关键的部分就是要把\(\prod_{j\neq i}(k-x_j)\)“展开”成一个关于\(k\)的求和的形式。我们先预处理出\(\prod_{j=0}^{n}(k-x_j)\)：把它写成一个以\(k\)为自变量的最高次数为\(n+1\)次的多项式。然后对于每个\(i\)，我们在预处理出的这个\(n+1\)次多项式的基础上，除以\((k-x_i)\)。
乘以\((k-x_i)\)和除以\((k-x_i)\)的操作都可以\(O(n)\)实现。具体来说，考虑两个生成函数：\(H(k)\)和\(G(k)\)。其中\(H(k)=G(k)(k-c)\)（请注意：\(c\)是在此处是一个常数，也就是前面式子里的\(x_i\)。而\(k\)是变量）。考虑生成函数的\(k^i\)项，可知：\(h_i=g_{i-1}-cg_i\)。那么，乘以\((k-c)\)相当于由\(G\)推\(H\)，除以\((k-c)\)相当于由\(H\)推\(G\)。根据\(h_i=g_{i-1}-cg_i\)直接递推即可。
总时间复杂度\(O(n^2)\)。
乘/除以\((k-x_i)\)的代码：
void mul(int f[],int c){//mul(k-c)
for(int i=n;i>=0;--i){
f[i]=(i?f[i-1]:0)-a*f[i];
}
}
void div(int f[],int c){//div(k-c)
int lst=0;
for(int i=n;i>=0;--i){
lst+=a*f[i+1];
swap(f[i],lst);
}
}
LOJ165: 动态加点，随时询问
朴素的做法中，加点是\(O(1)\)的；询问时按拉格朗日插值的公式模拟一遍，复杂度\(O(n^2)\)。因为有\(O(n)\)次询问，所以总时间复杂度\(O(n^3)\)。
发现朴素做法中，加点、查询的复杂度非常不平均：加点太快了、查询太慢了。于是我们想，能不能在加点的同时维护一些东西，使得两种操作的复杂度都变为\(O(n^2\log x)\)或\(O(n^2)\)。
考虑在询问时还是枚举\(i\)。我们预先维护好求积式中分母、分子的值。
- 分母：对于每个\(i\)，维护一个\(\text{fm}[i]=\frac{1}{\prod_{j\neq i}(x_i-x_j)}\)，这个可以在加点时顺便维护出来。
- 分子：在询问时，先特判\(k\)等于其中某个\(x_t\)的情况（直接输出\(y_t\)）。否则就求一个变量\(\text{fz}=\prod_{j=0}^{n}(k-x_j)\)，然后在枚举\(i\)时除以\((k-x_i)\)即可。
于是，一次询问的答案可以写为：\(\sum_{i=0}^{n}y_i\cdot \text{fm}[i]\cdot \text{fz}\cdot\frac{1}{k-x_i}\)。
如果在线回答询问，就需要\(O(\log x)\)求逆元。时间复杂度\(O(n^2\log x)\)。如果把询问离线下来，预处理所有逆元，则时间复杂度\(O(n^2)\)。两种做法的参考代码：在线(TLE 70分)，离线预处理逆元(AC 100分)。
\(x\)取值连续时
在\(x\)取值为连续的\(n+1\)个自然数时，我们可以把拉格朗日插值算法的时间复杂度优化为\(O(n)\)。
回到最开始的式子，即：\(f(k)=\sum_{i=0}^{n}y_i\prod_{j\neq i}\frac{k-x_j}{x_i-x_j}\)。
首先，我们可以预处理\(k-x_i\)的前缀积、后缀积。具体地，预处理出：\(\text{pre}[i]=\prod_{j=0}^{i}(k-x_j)\), \(\text{suf}[i]=\prod_{j=i}^{n}(k-x_j)\)。这样，上式中的分子就等于\(\prod_{j\neq i}(k-x_j)=\text{pre}[i-1]\cdot \text{suf}[i+1]\)。
然后考虑分母，即：\(\prod_{j\neq i}(x_i-x_j)\)。因为\(x\)取值连续，所以这个式子就等于\(\prod_{j\neq i}(i-j)\)。分\(j<i\)、\(j>i\)讨论，这个式子可以化为：\(i!\cdot(-1)^{n-i}\cdot(n-i)!\)。
综上所述，此时：
按这个式子直接模拟，时间复杂度\(O(n)\)。
自然数幂求和
给定\(n\), \(k\)，求：\(S(n,k)=\sum_{i=1}^{n}i^k\)。一般\(n\)很大（\(10^9\)级别）而\(k\)较小（\(10^5\)级别）。
其中一个经典的方法是把\(S(n,k)\)中\(n\)视为定值，然后写出其指数生成函数：\(\sum_{i=0}^{\infty}S(n,i)\frac{x^i}{i!}\)。再根据\(e^x=\sum_{i=0}^{\infty}\frac{x^i}{i!}\)进行一系列代换。最后用多项式求逆求出答案。时间复杂度\(O(k\log k)\)。这种方法的好处是可以一次性求出所有\(S(n,0\dots k)\)。
本文要介绍的拉格朗日插值法，一次只能求出一个\(S(n,k)\)，但是时间复杂度只有\(O(k)\)，而且更好写。
这种方法基于一个结论：
\(S(n,k)\)是以\(n\)为自变量的\(k+1\)次多项式。（请注意，这里与上一种方法恰恰相反，它把\(n\)看做变量、\(k\)看做定值）。
例如，\(\sum_{i=1}^{n}i^2=\frac{n(n+1)(2n+1)}{6}\)，这就是一个\(3\)次多项式。
因此我们可以带入\(k+2\)个点，然后用拉格朗日插值法直接求出\(S(n,k)\)。如果带入的\(k+2\)个点为\(1\dots k+2\)，则时间复杂度为\(O(k)\)。
至于预处理\(1^k,2^k,\dots ,(k+2)^k\)，可以用线性筛。在质数处暴力快速幂。因为质数只有\(O(\frac{k}{\log k})\)个。所以复杂度就是\(O(\frac{k}{\log k}\cdot \log k)=O(k)\)。
