拉格朗日插值

简介

对一个次数不超过 $n$ 的多项式，拉格朗日插值法可以在 $O(n^2)$ 的时间内，利用多项式的 $n+1$ 个取值，计算出多项式的各项系数。

多项式的系数表示法

众所周知一个 $k$ 次多项式可以被它的 $k+1$ 个系数唯一地确定，即：

$f(x)=\sum_{i=0}^ka_ix^i$

这种确定多项式的方法叫做多项式的**「系数表示法」**。

多项式的点值表示法

对于多项式 $f(x)$，如果 $f(x_0)=y_0$，不引起歧义的话，我们就可以说多项式 $f$ 在 $x_0$ 处的点值是 $y_0$。

还是众所周知地，一个 $k$ 次多项式同样可以被它的 $k+1$ 个点值所唯一地确定。

其中的道理不难讲：设这个多项式为

$f(x)=\sum_{i=0}^{k}a_ix^i$

那么，如果我们有了这个多项式上的 $k+1$ 个点 $(x_1,y_1),(x_2,y_2),\cdots,(x_{k+1},y_{k+1})$，就相当于有了一个 $k+1$ 次方程组：

$\begin{cases} a_0+a_1x_1+a_2x_1^2+\cdots+a_kx_1^k&=&y_1\\ a_0+a_1x_2+a_2x_2^2+\cdots+a_kx_2^k&=&y_2\\ \cdots&\cdots&\cdots\\ a_0+a_1x_{k+1}+a_2x_{k+1}^2+\cdots+a_kx_{k+1}^{k}&=&y_{k+1} \end{cases}$

这是一个由 $k+1$ 条 $k+1$ 元一次方程组成的一个方程组，并且不难证明在这 $k+1$ 条方程中，不存在两个等价的方程。

因此 $a_0,a_1,a_2,\cdots,a_{k}$ 这 $k+1$ 个未知数可以被唯一地解出来。也就是说，这 $k+1$ 个点同样唯一地确定了多项式

$f(x)=\sum_{i=0}^{k}a_ix^i$

这种确定多项式的方法叫做**「点值表示法」**。

系数表示法转点值表示法

直接 $O(n^2)$ 算。略。

对于一些特殊的点值，我们可以在更低的复杂度内算出这些点值。

点值表示法转系数表示法

首先由之前的推导过程，可以发先转换的过程实际上就是解一个 $k+1$ 元一次方程组。

因此可以利用高斯消元在 $O(k^3)$ 的时间内求解方程，确定其系数。

我们再来看 $O(k^2)$ 的拉格朗日插值。

首先有一个结论：

$f(x)\equiv f(a)\quad(\bmod(x-a))$

这是因为当 $x=a$ 时，有 $f(x)-f(a)=0$。

那么由我们初中就学过的因式定理，就有 $f(x)-f(a)\equiv0(\bmod(x-a))$，即 $f(x)\equiv f(a)\quad(\bmod(x-a))$。

分别将 $a$ 取 $x_1,x_2,\cdots,x_{k+1}$ ，那么 $f(a)$ 就取遍了 $y_1,y_2,\cdots,y_{k+1}$。

因此就得到一个关于 $f(x)$ 的方程组：

$\begin{cases} f(x)\equiv y_1(\bmod\ (x-x_1))\\ f(x)\equiv y_2(\bmod\ (x-x_2))\\ \cdots\\ f(x)\equiv y_{k+1}(\bmod\ (x-x_{k+1})) \end{cases}$

发现这个形式很像 中国剩余定理 对不对qwq

回忆一下中国剩余定理的过程：

1. 计算所有模数的积 $\text{Mul}=\prod\limits_{i=1}^{k+1}(x-x_i)$。

2. 对每个模数，计算 $m_i=\dfrac{\text{Mul}}{(x-x_i)}=\prod\limits_{i\neq j}(x-x_j)$。

3. 计算 $m_i$ 在模第 $i$ 个模数 $(x-x_i)$ 意义下的逆元：$m_i^{-1}=\prod\limits_{i\neq j}\dfrac{1}{x_i-x_j}$

4. 计算 $c_i=m_im_i^{-1}=\prod\limits_{i\neq j}\dfrac{x-x_j}{x_i-x_j}$

5. 线性同余方程组的解即为：$\sum\limits_{i=1}^{k+1}y_ic_i$，即：

$f(x)=\sum_{i=1}^{k+1}y_i\prod\limits_{i\neq j}\dfrac{x-x_j}{x_i-x_j}$

这就是 $f(x)$ 的表达式了qwq。

事实上，我们取一个点值 $(x_c,y_c)$，验证一下 $f(x_c)$，可以发现：

$f(x_c)=\sum_{i=1}^{k+1}y_i\prod_{i\neq j}\dfrac{x_c-x_j}{x_i-x_j}$

而当 $i\neq c$ 时，后面的那个 $\prod\limits_{i\neq j}\dfrac{x_c-x_j}{x_i-x_j}$ 中，总存在 $j$ 使得 $j=c$，那么 $x_c-x_j=0$，乘起来也是 $0$，那这一整项就都是 $0$。

这意味着那么大个求和符号中，实际上被计算的只有 $i=c$ 时的值，也就是： $y_c\prod\limits_{c\neq j}\dfrac{x_c-x_j}{x_c-x_j}=y_c$，因此确实有 $f(x_c)=y_c$。

于是 $f(x)$ 符合要求。我们只需要进行简单的二元多项式乘法，就可以计算出 $f(x)$ 的表达式。

不难发现这个表达式的计算复杂度是 $O(k^2)$ 的，那么我们就在 $O(k^2)$ 的时间内算出了 $f(x)$。

Coding 注意事项&技巧&实现

这里放上计算 $f(k)$ 的能过板子题的代码qwq。

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<iostream>

#define LL long long
const LL mod=998244353;

using namespace std;

LL n,k;
LL x[2005],y[2005];

LL ksm(LL a,LL b){
    LL ans=1ll,y=a;
    while(b){
        if(b&1)ans*=y,ans%=mod;
        y*=y;y%=mod;
        b>>=1;
    }
    return ans%mod;
}//快速幂，用于求逆元 

LL inv(LL x){
    return ksm(x,mod-2);
}//逆元qwq 

LL calc(LL t){
    LL ans=0;
    for(int i=1;i<=n;i++){
        LL p=y[i]%mod,q=1ll;
        for(int j=1;j<=n;j++){
            if(i==j)continue;
            p=p*(t-x[j])%mod;
            q=q*(x[i]-x[j])%mod;//直接照着公式算qwq
        }
        ans+=p*inv(q)%mod;
        ans%=mod;
    }
    return (ans%mod+mod)%mod;
}

int main(void){

    cin>>n>>k;
    for(int i=1;i<=n;i++)cin>>x[i]>>y[i];
    cout<<calc(k)<<endl;

    return 0;
}