数论 拉格朗日插值法(Lagrange Interpolation Polynomial)

· · 个人记录

零.简介

历史

在数值分析中，拉格朗日插值法是以法国十八世纪数学家约瑟夫·路易斯·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。拉格朗日插值法最早被英国数学家爱德华·华林于1779年发现，不久后（1783年）由莱昂哈德·欧拉再次发现。1795年，拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法，从此他的名字就和这个方法联系在一起。[1]

概念

一般地，若已知y=f(x)在乎不相同n+1个点x_0,x_1 \cdots ,x_n处的函数值y_0,y_1 \cdots ,y_n（即该函数通过(x_0,y_0),(x_1,y_1), \cdots , (x_n,y_n)这n+1个点），则可以考虑构造一个次数不超过n的多项式y=Pn(x)使其满足:

Pn(x_i)=y_i (i \in (0,1,\cdots n))

任意一个实数\sigma(\sigma \notin (x_0 , x_1 \cdots x_n))都可以用Pn(\sigma)的值作为f(\sigma)的准确值的近似值，这个方法叫作拉格朗日插值法

我们将上面Pn的表达式称为插值条件（准则）

定理: 满足插值条件的、次数不超过n的多项式是存在而且是唯一的 [1]

用人话转述

如果我们知道一条函数f(x)过n+1过不同的点

那么我们就可以唯一地确定一条次数不超过n的函数g(x)来作为f(x)的近似函数

并且可以用g(x_0)来作为f(x_0)的近似值 其中x_0为任意实数

而这个方法就被称为拉格朗日插值法

一.拉格朗日插值法

1.思路

我们令D=\{ 0,1,\cdots ,n\}

在拉格朗日插值法中我们把最终得到的函数g(x)分为n+1个部分

我们把这些部分记作h_i(i\in D)

试想如果我们让x=x_i时可以使h_i=y_i而其他h都是0

那对应的函数是不是就是我们要求的函数了呢?

从x_i的角度来看是如此

那我们从h_i的角度来看就是:

(1).x =x_j(j\in D \& j\not= i)$时$h_i=0 (2).x=x_i$时$h_i=y_i

先不妨假设n=4且i=2来推

根据(1)可以得到h_2=k\times (x-x_0)(x-x_1)(x-x_3)

再把(2)代入:k\times(x_2-x_0)(x_2-x_1)(x_2-x_3)=y_2

所以k=\dfrac{y_2}{(x_2-x_0)(x_2-x_1)(x_2-x_3)}

所以h_2=\dfrac{(x-x_0)(x-x_1)(x-x_3)}{(x_2-x_0)(x_2-x_1)(x_2-x_3)}\times y_2

推广到一般的情况就是:

\because (1) \therefore h_i=k\times \prod_{j\in D \& j\not= i} (x-x_j) \because (2)\therefore k \times \prod_{j\in D \& j\not= i}(x_i-x_j)=y_i k=\dfrac{y_i}{\prod(x_i-x_j)}(j\in D \& j\not= i) \therefore h_i=\dfrac{\prod (x-x_j)}{\prod(x_i-x_j)}\times y_i (j\in D \& j\not= i)

那么新得到的函数g(x)就是

g(x)=\sum_{i=0}^n y_i\prod_{j \in D \& j \not=i}\dfrac{(x-x_j)}{x_i-x_j}

在这里举一个例子 n=3时得到的函数就是:

g(x)=\dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}\times y_0+\dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} \times y_1+\dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}\times y_2

代x_0,x_1,x_2进去看一下是不是满足?

2.例题\&相当于模板

【模板】拉格朗日插值

【题目描述】

由小学知识可知 n 个点 (x_i,y_i)可以唯一地确定一个多项式y=f(x)。

现在，给定这 n 个点，请你确定这个多项式，并求出 f(k) \bmod 998244353的值。

【输入格式】

第一行两个整数n(n≤ 2000),k

接下来n行，第i行两个整数x_i,y_i

【输出格式】

一行一个整数，表示f(k) \bmod 998244353的答案

样例输入

3 100
1 4
2 9
3 16

样例输出

10201

正如题目所言是一道模板题

我们就O(n^2)枚举i和j求答案即可

code

#include<cstdio>
#include<cctype>
const int mod=998244353;
inline int R(){
    int r=0;char c=getchar();
    while(!isdigit(c)) c=getchar();
    while(isdigit(c)) r=(r<<1)+(r<<3)+(c^48),c=getchar();
    return r;
}
int n,k;
long long x[2005],y[2005];
inline long long inv(long long a){//求逆元
    int k=mod-2;
    long long node=a,ans=1;
    while(k){
        if(k&1) ans=(node*ans)%mod;
        node=(node*node)%mod;
        k>>=1;
    }
    return ans;
}
int main(){
    n=R();k=R();
    for(register int i=1;i<=n;i++){
        x[i]=R();y[i]=R();
    }
    long long ans=0;
    for(register int i=1;i<=n;i++){
        long long sum=y[i];//sum即是h_i
        for(register int j=1;j<=n;j++){
            if(j==i) continue;
            sum=(sum*(k-x[j]))%mod;
            sum=(sum*inv(x[i]-x[j]))%mod;
        }
        ans+=sum;ans=(ans+mod)%mod;
    }
    printf("%lld\n",ans);
    return 0;
}

二.重心拉格朗日插值

g(k)=\sum_{i=0}^n y_i \prod_{j \in D \& j \not=i}\dfrac{(k-x_j)}{x_i-x_j}

再看一遍这个式子就发现

在程序中我们重复计算了很多次k-x_j

我们就新设一个参数g=\prod (k-x_i)

我们就把式子中的(k-x_j)提出来 但是还剩一个k-x_i

g(k)=g\sum_{i=0}^n( \dfrac{y_i}{k-x_i} \prod_{j \in D \& j \not=i }(\dfrac{1}{x_i-x_j}))

再把只和i有关的提出来，再设t_i=\dfrac{y_i}{\prod_{j\in D \& j \not=i} (x_i-x_j)}

于是式子就变成了:

g(k)=g \sum_{i=0}^n \dfrac{t_i}{k-x_i}

看下这个式子是不是简洁多了

于是我们每次只求出t_i就可以了

在这种情况下上面的模板题代码就是:

#include<cstdio>
#include<cctype>
const int mod=998244353;
inline int R(){
    int r=0;char c=getchar();
    while(!isdigit(c)) c=getchar();
    while(isdigit(c)) r=(r<<1)+(r<<3)+(c^48),c=getchar();
    return r;
}
inline long long fastm(int a,int k){
    long long ans=1,node=a;
    while(k){
        if(k&1) ans=(ans*node)%mod;
        node=(node*node)%mod;
        k>>=1;
    }
    return ans;
}
int n,k;
int x[2005],y[2005];
long long g=1,t[2005],ans;
int main(){
    n=R();k=R();
    for(register int i=1;i<=n;i++){
        x[i]=R();y[i]=R();
        g=g*(k-x[i])%mod;
    }
    for(register int i=1;i<=n;i++){
        t[i]=y[i];
        for(register int j=1;j<=n;j++){
            if(i==j) continue;
            t[i]=(t[i]*fastm(x[i]-x[j],mod-2))%mod;
        }
    }
    for(register int i=1;i<=n;i++){
        ans+=(t[i]*fastm(k-x[i],mod-2))%mod;
        ans%=mod;
    }
    printf("%lld\n",(ans*g%mod+mod)%mod);
    return 0;
}

三.总结

拉格朗日插值法真的是一个很神奇的东西

这个方法启示我们没有最难的问题，只有更巧妙的构造

我们可以通过它找出题目中涉及到未知的高阶函数

当然我们很多时候我们也没办法确定一个问题是否是一条高阶函数

更多时候会在式子题会有应用

Finally，谢谢大家

推荐习题

参考资料

[1] .拉格朗日插值法[OL].百度百科.2020-08-21/2021-4-30.