拉格朗日插值

· · 算法·理论

以前听到拉插感觉很高深，很不可做，但是真正自己做起来发现也没想象中的那么困难喵。

起因是 NOIP2020 的 T4 要拉插，所以跑去学的。

P4781 【模板】拉格朗日插值

插值

插值，指的是通过已知的一些点去推算另一些点的值，简单来说就是通过一些点去拟合出一个函数来。一般我们碰到的都是多项式插值，即用多项式去拟合一个函数。

引理：用 n+1 个点能唯一拟合出一个 n 阶多项式。

感性理解即可（？）

Lagrange 插值法

对于多项式插值来说，就是要构造一个 n+1 阶函数 f(x) 满足过点 P_1(x_1,y_1),P_2(x_2,y_2)\dots P_n(x_n,y_n)。我们考虑构造 n 个函数 g_1(x),g_2(x)\dots g_n(x)，满足：

g_i(x)=\begin{cases} P_j(x_j,0) & j\neq i \\P_i(x_i,y_i) & j=i \end{cases}

简单来说就是对于每个函数来讲，保证过第 i 个点，而且对于其他的 x 来讲都要取到 0。这样的话就不难得出最终的 f(x)=\sum\limits_{i=1}^n g_i(x)。

接下来考虑如何构造函数 g_i。由于对于 x=x_j(j\neq i) 都取 0，那么不难发现可以用一个系数 \prod_{j\neq i}\limits x-x_j 来表示。然后考虑如何使得在 P_i 的时候能取到对应的点。很简单的就是可以乘上一个 \frac{y_i}{\prod\limits_{j\neq i} x_i-x_j}。这样就能使得这个函数满足上面的条件了。总的来说，g_i(x)=y_i\times \prod\limits_{j\neq i}\frac{x-x_j}{x_i-x_j}。所以说 Lagrange 插值的形式就是：

f(x)=\sum\limits_{i=1}^n y_i\times \prod\limits_{j\neq i}\frac{x-x_j}{x_i-x_j}

不难吧（）时间复杂度 O(n^2)。

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define ull unsigned long long
#define pii pair<int,int>
#define ls (root << 1)
#define rs (root << 1 | 1)
#define mid ((l + r) >> 1)
const int MAX = 2e5 + 5;
const int mod = 998244353;
const int INF = 0x7fffffff;
inline int read(){
    int x=0,y=1;char c=getchar();
    for(;c<'0'||c>'9';c=getchar())if(c=='-')y=-1;
    for(;c>='0'&&c<='9';c=getchar())x=(x<<1)+(x<<3)+(c^48);
    return x*y;
}
inline int qpow(int a,int b){int res=1;for(;b;b>>=1,a=a*a%mod)if(b&1)res=res*a%mod;return res;}
int n,k,x[MAX],y[MAX],ans;
signed main()
{
//  freopen(".in","r",stdin),freopen(".out","w",stdout);
    n=read(),k=read();
    for(int i=1;i<=n;i++) x[i]=read(),y[i]=read();
    for(int i=1;i<=n;i++)
    {
        int pro1=1,pro2=1;
        for(int j=1;j<=n;j++) if(i!=j) pro1=pro1*(k-x[j])%mod,pro2=pro2*(x[i]-x[j])%mod;
        ans=(ans+y[i]*pro1%mod*qpow(pro2,mod-2))%mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}

例题：CF622F

F. The Sum of the k-th Powers

题目大意：用 O(k\log k) 的复杂度求解 \sum\limits_{i=1}^n i^k。

想不到吧，这题要用拉插喵）

这个柿子可以理解为一个 n+1 次的多项式。可以感性理解，也可以这么理解：k=1 时是一个等差数列，化简之后就是一个两阶的多项式；k\ge 2 的时候也可以差不多这样递推理解（？）

这样知道之后就好拉插了。求出 前 k+2 项，然后就直接拉插就好了（吗），显然不是。由于正常的拉插的复杂度是 O(k^2) 级别的，肯定过不了啊。但是根据这个题，注意到它的 x 的值是连续的，这启发我们似乎可以用一些前缀、预处理的手段去实现单次计算 O(1)，然后乘上求逆元、快速幂的 O(\log k) 的复杂度就可以达到 O(k\log k) 了。

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define ull unsigned long long
#define pii pair<int,int>
#define ls (root << 1)
#define rs (root << 1 | 1)
#define mid ((l + r) >> 1)
const int MAX = 1e6 + 10;
const int mod = 1e9+7;
const int INF = 0x7fffffff;
inline int read(){
    int x=0,y=1;char c=getchar();
    for(;c<'0'||c>'9';c=getchar())if(c=='-')y=-1;
    for(;c>='0'&&c<='9';c=getchar())x=(x<<1)+(x<<3)+(c^48);
    return x*y;
}
int n,k,fac[MAX],pro=1;
inline int qpow(int a,int b){int res=1;for(;b;b>>=1,a=a*a%mod)if(b&1)res=res*a%mod;return res;}
signed main()
{
//  freopen(".in","r",stdin),freopen(".out","w",stdout);
    n=read(),k=read();
    int ans=0,y=0;
    if(n<=k+2) {for(int i=1;i<=n;i++) ans=(ans+qpow(i,k))%mod;printf("%lld\n",ans);return 0;}
    fac[0]=1;
    for(int i=1;i<=k+2;i++) fac[i]=fac[i-1]*i%mod;
    for(int i=1;i<=k+2;i++) pro=pro*(n-i)%mod;
    for(int i=1;i<=k+2;i++) y=(y+qpow(i,k))%mod,
        ans=(ans+y*pro%mod*qpow((n-i)*fac[i-1]%mod*fac[k+2-i]%mod,mod-2)%mod*((k+2-i)%2?-1:1))%mod;
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}

oi-wiki 传送门