引入

在初中时，我们都学过一次函数，一定也都做过这种题：

给定两个点，求经过这两个点的一次函数的解析式。

我们只需要待定系数法，然后把点坐标带进去解方程即可。

那么推广一下，可以得到结论：

在平面直角坐标系中，$n+1$ 个点能唯一确定一个 $n$ 次多项式。

但如果我们暴力高斯消元，复杂度是 $O(n^3)$ 的。

而拉格朗日插值法可以在 $O(n^2)$ 的复杂度内解决这个问题。

单点求值

设 $f(x)$ 为待求多项式，$(x_i,y_i)$ 为第 $i$ 个点的坐标。

以下均保证 $x_i$ 互不相同。

设拉格朗日基本多项式为

$$g_j(x)=\prod_{i=0,i\neq j}^n \frac {x-x_i}{x_j-x_i}$$

这个多项式十分巧妙。

可以发现，$\forall \ i\in[0,n],g_j(x_i)=[i=j]$。

根据基本多项式，我们可以构造出 $f(x)$：

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

根据基本多项式的性质，我们可以得到 $f(x_i)=y_i$，也就是恰好经过这 $n+1$ 个点。

那么单点求值时，只需要把那个点带进去即可。

复杂度 $O(n^2)$

事实上应该线性求逆元才是严格 $O(n^2)$ 的，但我懒了，所以写的是 $O(n^2\log n)$ 的。

例题：【模板】拉格朗日插值

#include"bits/stdc++.h"
#define re register
#define int long long
using namespace std;
const int maxn=2e3+10,mod=998244353;
int n,k,ans;
int x[maxn],y[maxn];
inline int qpow(int a,int b){
    int res=1;
    while(b){
        if(b&1) res=res*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return res;
}
inline int Inv(int x){return qpow(x,mod-2);}
signed main(){
#ifndef ONLINE_JUDGE
    freopen("1.in","r",stdin);
    freopen("1.out","w",stdout);
#endif
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>k;
    for(re int i=1;i<=n;++i) cin>>x[i]>>y[i];
    for(re int i=1;i<=n;++i){
        int s1=y[i],s2=1;
        for(re int j=1;j<=n;++j){
            if(i!=j){
                s1=s1*((k-x[j])%mod+mod)%mod;
                s2=s2*((x[i]-x[j])%mod+mod)%mod;
            }
        }
        ans=(ans+s1*Inv(s2)%mod)%mod;
    }
    cout<<ans;
    return 0;
}

求出系数

现在我们已经可以做到 $O(n^2)$ 复杂度的单点求值了。

但这还不够，因为我想具体求出来这个多项式。

把上面那个东西拿过来：

$$f(x)=\sum_{i=0}^n y_i\prod_{j=0,j\neq i}^n\frac {x-x_j}{x_i-x_j}$$

提出常数部分：

$$a_i=\frac{y_i}{\prod_{j=0,j\neq i}^n(x_i-x_j)}$$

然后设 $g(k)=\prod_{i=0}^n(k-x_i)$，则

$$f(x)=\sum_{i=0}^n a_i\frac {g(x)}{x-x_i}$$

我们可以暴力的把 $g(x)$ 这个多项式算出来，这是 $O(n^2)$ 的。

然后对于每一项，我们暴力模拟长除法，总复杂度是 $O(n^2)$ 的。

综上，我们在 $O(n^2)$ 的复杂度内得到了多项式的系数。

#include"bits/stdc++.h"
#define re register
#define int long long
using namespace std;
const int maxn=2e3+10,mod=998244353;
int n,k;
int x[maxn],y[maxn];
int f[maxn],g[maxn];
int a[maxn];
int tmp[maxn];
inline int qpow(int a,int b){
    int res=1;
    while(b){
        if(b&1) res=res*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return res;
}
inline int Inv(int x){return qpow(x,mod-2);}
inline void pre(){
    for(re int i=1;i<=n;++i){
        int res=1;
        for(re int j=1;j<=n;++j){
            if(i!=j) res=res*(x[i]-x[j]+mod)%mod; 
        }
        a[i]=y[i]*Inv(res)%mod;
    }
    g[0]=1;
    for(re int i=1;i<=n;++i,swap(g,tmp)){
        tmp[0]=0;
        for(re int j=1;j<=i;++j) tmp[j]=g[j-1];
        for(re int j=0;j<=i;++j) tmp[j]=(tmp[j]+g[j]*(-x[i]+mod)%mod)%mod;
    }
}
inline void Lagrange(){
    for(re int i=1;i<=n;++i){
        int lst=0,inv=Inv(-x[i]+mod);
        for(re int j=0;j<n;++j){
            tmp[j]=(g[j]-lst+mod)*inv%mod;
            f[j]=(f[j]+a[i]*tmp[j]%mod)%mod,lst=tmp[j];
        }
    }
}
inline int calc(int x){
    int ans=0,res=1;
    for(re int i=0;i<n;++i){
        ans=(ans+res*f[i]%mod)%mod;
        res=res*x%mod;
    }
    return ans;
}
signed main(){
#ifndef ONLINE_JUDGE
    freopen("1.in","r",stdin);
    freopen("1.out","w",stdout);
#endif
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>k;
    for(re int i=1;i<=n;++i) cin>>x[i]>>y[i];
    pre();
    Lagrange();
    cout<<calc(k);
    return 0;
}

x 取值连续时的优化

此时式子变为

$$f(k)=\sum_{i=0}^n y_i\prod_{j=0,i\neq j}^n \frac{k-i}{i-j}$$

对于分子，我们维护前缀积和后缀积。

对于分母，可以发现这是一个阶乘形式。

所以式子变为

$$f(k)=\sum_{i=0}^n y_i\frac {pre_{i-1}\times suf_{i+1}}{(i-1)!\times (n-i)!}$$

当 $n-i$ 为奇数时，分母应为负号。

复杂度 $O(n)$

例题：求 $\sum_{i=1}^n i^k$。

The Sum of the k-th Powers

可以证明，$\sum_{i=1}^n i^k$ 为 $k+1$ 次多项式。

所以直接取 $k+2$ 个点把它插出来即可。

取值连续，所以可以做到 $O(n)$。

#include"bits/stdc++.h"
#define re register
#define int long long
using namespace std;
const int maxn=1e6+10,mod=1e9+7;
int n,m,k,ans;
int fac[maxn],pre[maxn],suf[maxn];
inline int qpow(int a,int b){
    int res=1;
    while(b){
        if(b&1) res=res*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return res;
}
inline int Inv(int x){return qpow(x,mod-2);}
signed main(){
#ifndef ONLINE_JUDGE
    freopen("1.in","r",stdin);
    freopen("1.out","w",stdout);
#endif
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>k;m=k+2;
    pre[0]=fac[0]=suf[m+1]=1;
    for(re int i=1;i<=m;++i) pre[i]=pre[i-1]*(n-i)%mod;
    for(re int i=m;i>=1;--i) suf[i]=suf[i+1]*(n-i)%mod;
    for(re int i=1;i<=m;++i) fac[i]=fac[i-1]*i%mod;
    int y=0;
    for(re int i=1;i<=m;++i){
        y=(y+qpow(i,k)%mod)%mod;
        int s1=pre[i-1]*suf[i+1]%mod;
        int s2=((m-i)&1?-1ll:1ll)*fac[i-1]*fac[m-i]%mod;
        ans=(ans+(y*s1%mod*Inv(s2)%mod)%mod+mod)%mod;
    }
    cout<<ans;
    return 0;
}