洛谷P4781 【模板】拉格朗日插值

https://www.luogu.org/problemnew/show/P4781

https://www.cnblogs.com/zj75211/p/8028165.html讲的很清楚了

（这里有隐藏的转载&备份）

今天做题需要用到插值法, 就简单入门了一下, 同时记了一点浅显的东西于此。

定理

对于给定的 N+1个 (x坐标不相等的二维平面上的)点，存在唯一一个最高项次数不超过 N 的多项式 Y=P(X)其图像经过这N+1个点。

作用

求出的插值函数P(X)用于估计原函数F[X]。

(但如果原函数可以由多项式表示(既不是一个相关函数),且告诉了图像上其最高次项次数+1个点，就可以通过插值法得出该函数F[x],即求得的插值函数P[x]=F[x])

 求法（对于给定的一组N+1个点(X₀,y₀),(X₁,y₁)......(X_N,Y_N)）

1.直接法

设出N次方程，得到N+1个线性方程，求解出多项式的系数即可。

（高斯消元哦，N³哦，很慢哦）

2.拉格朗日插值法(可以叫它构造法么？)

先构造出一组(N+1个)基函数 l₀(X) , l₁(X) , ...... , l_N(X)

使得 l_i (X)只经过(X_i , 1)，且在X=X_j ( j ≠ i )时，l_i (X_j)=0。

l_i(X)的构造如上。

显然，上式在 X=X_j ( j ≠ i )时，l_i ( X_j )=0

同时在 X=X_i 时，l_i(X_i)=1，即 y_i l_i ( X_i )=y_i

然后插值函数 P(X)即为 y_i l_i ( X) 的线性相加：

即 P(X)=y₀ l₀ ( X)+y₁ l₁ ( X)+y₂ l₂ ( X)+……+y_N l_N ( X)

以下四个点为例：

P(X)即为我们所求的函数。

复杂度O(N²)

想要了解更多的话就看看这里：★

---

拉格朗日插值

设给定点为$(x_0,y_0)$到$(x_{n-1},y_{n-1})$，求经过它们的多项式F

公式：$F(x)=sum_{i=0}^{n-1}y_iprod_{j=0,j e i}^{n-1}frac{x-x_j}{x_i-x_j}$

对于此题，不需要真的把多项式求出来，直接在计算过程中把给定数带入求值就行（$O(n^2)$）

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
typedef long long ll;
typedef unsigned long long ull;
const int md=998244353;
int n,K,ans;
int x[2011],y[2011];
int poww(int a,int b)
{
    int ans=1;
    for(;b;b>>=1,a=ll(a)*a%md)
        if(b&1)
            ans=ll(a)*ans%md;
    return ans;
}
int main()
{
    int i,j,t1,t2;
    scanf("%d%d",&n,&K);
    for(i=0;i<n;++i)
        scanf("%d%d",x+i,y+i);
    for(i=0;i<n;++i)
    {
        t1=t2=1;
        for(j=0;j<n;++j)
            if(j!=i)
            {
                t1=ll(t1)*(K-x[j]+md)%md;
                t2=ll(t2)*(x[i]-x[j]+md)%md;
            }
        ans=(ans+ll(t1)*poww(t2,md-2)%md*y[i])%md;
    }
    printf("%d
",ans);
    return 0;
}

真的有必要求出多项式的话，可以先算出全部$(x-x_i)$的积，每次把它除以一项，来得到各个$l_i(x)$

这个“除以一项”怎么搞？

设除完$(x-x_i)$的结果为A(x)，各项系数为$a_i$，除之前为B(x)，各项系数为$b_i$，令$y=x_i$

则$b_{n+1}=a_n$，$b_n=a_{n-1}-a_ny$，$b_{n-1}=a_{n-2}-a_{n-1}y$，...

$a_n=b_{n+1}$，$a_{n-1}=b_n+a_ny$，$a_{n-2}=b_{n-1}+a_{n-1}y$，...

还是$O(n^2)$

当然，此题这么做要慢一点；不过这个在某些情况下有用，由于多项式求值可以O(n)完成

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
using namespace std;
#define fi first
#define se second
#define mp make_pair
#define pb push_back
typedef long long ll;
typedef unsigned long long ull;
const int md=998244353;
int n,K,ans;
int x[2011],y[2011];
int an1[2011],tt1[2011],tt2[2011];
//tt1维护所有(x-x_i)的积，an1维护答案
int poww(int a,int b)
{
    int ans=1;
    for(;b;b>>=1,a=ll(a)*a%md)
        if(b&1)
            ans=ll(a)*ans%md;
    return ans;
}
int main()
{
    int i,j,t2;
    scanf("%d%d",&n,&K);
    for(i=0;i<n;++i)
        scanf("%d%d",x+i,y+i);
    tt1[0]=1;
    for(i=0;i<n;++i)
    {
        tt1[i+1]=tt1[i];
        for(j=i;j>=1;--j)
            tt1[j]=(ll(tt1[j])*(md-x[i])+tt1[j-1])%md;
        tt1[0]=ll(tt1[0])*(md-x[i])%md;
    }
    for(i=0;i<n;++i)
    {
        tt2[n-1]=tt1[n];
        for(j=n-2;j>=0;--j)
            tt2[j]=(tt1[j+1]+ll(tt2[j+1])*x[i])%md;
        t2=1;
        for(j=0;j<n;++j)
            if(j!=i)
                t2=ll(t2)*(x[i]-x[j]+md)%md;
        t2=ll(poww(t2,md-2))*y[i]%md;
        for(j=0;j<n;++j)
            an1[j]=(an1[j]+ll(t2)*tt2[j])%md;
    }
    ans=0;
    for(j=n-1;j>=0;--j)
        ans=(ll(ans)*K+an1[j])%md;
    printf("%d
",ans);
    return 0;
}

---

重心拉格朗日插值

https://blog.csdn.net/qq_35649707/article/details/78018944

https://zh.wikipedia.org/wiki/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E6%8F%92%E5%80%BC%E6%B3%95

设$l(x) = prod_j(x-x_j)$

$omega_i= frac{1}{prodlimits_{j ot = i}(x_i-x_j)}$

则$F(x)=l(x)sum_ifrac{omega _i}{x-x_i}y_i$（的确是对的，过程省略了）

这个东西的好处是，在增加一个新点的坐标时，只需要O(n)的时间来更新解，不需要$O(n^2)$重新计算整个解（具体方法先鸽着）

还有一个公式（重心拉格朗日插值（第二型））

将$g(x)=1$在对原多项式给出的所有x坐标处用以上方法插值，得到$l(x)sum_ifrac{omega _i}{x-x_i}=1$（并不知道怎么想出来的）

所以$F(x)=frac{F(x)}{1}=frac{sum_ifrac{omega _i}{x-x_i}y_i}{sum_ifrac{omega _i}{x-x_i}}$（暂时不知道有什么特别的用处）

其他资料https://www.cnblogs.com/ECJTUACM-873284962/p/6833391.html