【ACM数论】拉格朗日插值法

拉格朗日插值法

1.拉格朗日插值法的由来

​ 拉格朗日插值法究竟解决的是什么问题喃？

问题：在二维图像上给出$n$个点，如何得到一个通过这$n$个点的一条线？

​ 例如：已知以下几个点$(0,1)、(1,0)、(2,1)$,问这条曲线是什么？

显然易见对于这三个点，是可以用多项式画出这三个点的，因此我们可以假设这个二次多项式为
$$y=a_0+a_{1}x+a_2{x}^2$$
我们带入这三个点$(x_1,y_1),(x_2,y_2),(x_3,y_3)$

可以得到以下方程组
$$\{\begin{array}{l}{y}_1=a_0+a_1{x}_1+a_2{x}_1^2\\ y_2=a_0+a_1{x}_2+a_2{x}_2^2\\ y_3=a_0+a_1{x}_3+a_2{x}_3^2\end{array}$$
对于这样一个方程组，我们可以用高斯消元解出系数。

但是拉格朗日是如何思考这个问题的喃？

拉格朗日的思考：

拉格朗日从另外的角度思考了这个问题，在他看来可以通过三根二次曲线相加来达到目标。那这是怎么的三根二次曲线呢？

这三根曲线只需满足一些条件即可：
$$\begin{gathered} 曲线1：y_1{f}_1(x)可以保证，在x_1点处，取值为y_1,其余两点取值为0 \\ 曲线2：y_2{f}_1(x)可以保证，在x_2点处，取值为y_2,其余两点取值为0 \\ 曲线3：y_3{f}_1(x)可以保证，在x_3点处，取值为y_3,其余两点取值为0 \end{gathered}$$
则所求曲线为：
$$f(x)=y_1{f}_1(x)+y_2{f}_2(x)+y_3{f}_3(x)$$
而这样一条曲线正好可以过给的三个点

2、拉格朗日插值法的推导

对于过每一个点$(x_i,y_i)$的曲线的函数设为$y_j{f}_i(x_j)$

因此需要满足的条件为：
$$f_i(x_j)=\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array}$$
因此我们最后构造出$f_i(x_j)$即可

显而易见，以下这个函数满足条件
$$f_i(x)=\frac{(x-x_1)(x-x_2)...(x-x_n)}{(x_i-x_1)(x_i-x_2)...(x_i-x_n)}(i\ne j)$$
因此
$$f(x)=\sum _{i=1}^n{y}_i{f}_i(x)=\sum _{i=1}^n{y}_i{\Pi }_{i\ne j}\frac{x-x_j}{x_i-x_j}$$
上模板题：https://www.luogu.com.cn/problem/P4781

#include<iostream>
#include<algorithm>
#include<cmath>
#include<string.h>
using namespace std;
#define int long long
const int maxn=1e4+10;
int x[maxn];
int y[maxn];
const int mod =998244353;
int qp(int x,int n){
    int ans=1;
    while(n){
        if(n&1){
            ans=ans*x%mod;
        }
        x=x*x%mod;
        n>>=1;
    }
    return ans;
}
int inv(int x){
    return qp(x,mod-2);
}
int lagrange(int num,int n)
{
    int ans=0;
    for(int k=1;k<=n;k++){
        int s1=y[k];
        int s2=1;
        for(int i=1;i<=n;i++){
            if(i==k)continue;
            s1=s1*(num-x[i])%mod;
            s2=s2*(x[k]-x[i])%mod;
        }
        ans=(ans+s1*inv(s2))%mod;
    }
    ans=(ans%mod+mod)%mod;
    return ans;
}

signed main()
{
    int n,k;
    cin>>n>>k;
    for(int i=1;i<=n;i++){
        cin>>x[i]>>y[i];
    }
    cout<<lagrange(k,n)<<endl;
    return 0;
}

当然这样的时间复杂度$O(n^2)$

在有些情况下不能做到最优，因此对于一些情况我们可以进行优化

3、横坐标为连续的自然数（1,2,3…）时，可以进行优化

由于横坐标为连续的自然数:

故插值公式可以优化为：
$$f(x)=\sum _{i=1}^{n+1}{y}_i\frac{{\Pi }_{j=1}^{n+1}(x-j)}{(x-j)(-1{)}^{n+1-i}(i-1)!(n+1-i)!}$$
这样一来，可以通过预处理后缀积、前缀积然后代入这个式子，所以这种情况时间复杂度为$O(n)$

例题：求自然数的k次幂的和
$$\begin{gathered} 求\sum _{i=1}^n{i}^k \\ 思路：首先我们通过这个式子可以判断出这个式子最终一定是一个k+1次的多项式 \\ 因此我们只需要k+2个点 \\ 那么接下来如何设计这个函数喃？我们可以把所求看为f(n)=\sum _{i=1}^n{i}^k \\ 那么每一个点为(j,\sum _{i=1}^j{i}^k) \\ 因此我们预处理出k+2个点，然后将这k+2个点代入插值公式即可 \end{gathered}$$
相同题意的例题：https://codeforces.com/problemset/problem/622/F

代码：

#include <iostream>
#include <algorithm>
#include <cmath>
#include <string.h>
#define int long long
using namespace std;
const int maxn = 1e6 + 10;
const int mod = 1e9 + 7;
int x[maxn];
int y[maxn];
int qp(int a, int b)
{
    int ans = 1;
    while (b)
    {
        if (b & 1)
        {
            ans = ans * a % mod;
        }
        a = a * a % mod;
        b >>= 1;
    }
    return ans;
}
int inv(int x)
{
    return qp(x, mod - 2);
}
int fac[maxn];
void init_fac()
{
    fac[0] = 1;
    fac[1] = 1;
    for (int i = 2; i < maxn - 5; i++)
    {
        fac[i] = fac[i - 1] * i % mod;
    }
}
int back;
int lagrange(int n, int k)
{
    int ans = 0;
    for (int i = 1; i <= k + 2; i++)
    {
        int s1 = back * y[i] % mod;
        int s2 = ((k + 2 - i) % 2 == 0 ? 1 : -1) * (n - i);
        s2 = s2 * fac[i - 1] % mod;
        s2 = s2 * fac[k + 2 - i] % mod;
        ans = (ans + s1 * inv(s2) % mod) % mod;
    }
    ans = (ans % mod + mod) % mod;
    return ans;
}
signed main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    int n, k;
    cin >> n >> k;
    init_fac();
    for (int i = 1; i <= k + 2; i++)
    {
        x[i] = i;
        y[i] = (y[i - 1] + qp(i, k)) % mod;
    }
    if (n <= k+2)
    {
        cout << y[n] << endl;
    }
    else
    {
        back = 1;
        for (int j = 1; j <= k + 2; j++)
        {
            back = back * (n - j) % mod;
        }
        cout << lagrange(n, k) << endl;
    }
    return 0;
}