[学习笔记]拉格朗日插值

对于拉格朗日插值的了解始于知乎上的数列问题：

问： \(1,3,5,7\) 的下一项是什么啊？
答：根据拉格朗日插值公式，可以显然地构造函数

\[\large f(x)=\dfrac{18111}{2}x^4-90555x^3+\dfrac{633885}{2}x^2-452773x+217331 \]

使得 \(f(1)=1,f(2)=3,f(3)=5,f(4)=7\) ，不难发现 \(f(5)=114514\)
也就是说原数列下一项为 \(114514\)

好了说正经的
众所周知，\(n+1\) 个 \(x\) 坐标不同的点可以确定唯一的最高为 \(n\) 次的多项式。
现在给你 \(n+1\) 个点，让你求这个多项式在 \(k\) 处的值
不妨设 \(n\) 次多项式为 \(f(x)\) ，给定 \(n\) 个点的坐标为 \((x_i,y_i)\)
我们有拉格朗日插值公式

\[\Large f(x)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j=1,j\neq i}^{n+1}\frac{x-x_j}{x_i-x_j} \]

首先这个是显然正确的：将任意 \(x_i\) 代入，除了第 \(i\) 项的其它所有东西都会因为分子有 \(-x_i\) 而被消去
没了

P4781 【模板】拉格朗日插值

点击查看代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<string>
#include<iostream>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099511627776,mod=998244353;
int n,k;
int x[WR],y[WR];
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<3)+(s<<1)+ch-'0';
        ch=getchar();
    }
    return s*w;
}
int quick_pow(int a,int b){
    int base=a,res=1;
    while(b){
        if(b&1) res=res*base%mod;
        base=base*base%mod;
        b>>=1;
    }
    return res;
}
int Lagrange(){
    int ans=0;
    for(int i=1;i<=n;i++){
        int up=1,dn=1;
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            up=up*(k-x[j])%mod;
            dn=dn*(x[i]-x[j])%mod;
        }
        ans=(ans+y[i]*up%mod*quick_pow(dn,mod-2)%mod)%mod;
    }
    while(ans<0) ans+=mod;
    return ans;
}
signed main(){
    n=read(),k=read();
    for(int i=1;i<=n;i++) x[i]=read(),y[i]=read();
    printf("%lld\n",Lagrange());
    return 0;
}

如果已知点的横坐标是连续整数，我们可以做到 \(\mathcal{O}(n)\) 插值。
此时 \(x_1=1~,~x_2=2~,~\cdots~,~x_n=n\)
注意到

\[\begin{aligned} f(x)&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{i\neq j}\frac{x-x_j}{x_i-x_j}\\ &=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{i\neq j}\frac{x-j}{i-j}\\ \end{aligned} \]

那么显然地，分子为

\[\frac{\prod\limits_{k=1}^{n+1}(x-k)}{x-i} \]

分母可以拆成两段算，为

\[(-1)^{n-i+1}(i-1)!(n+1-i)! \]

至此，原式化为

\[f(x)=\sum\limits_{i=1}^{n+1}y_i\frac{\prod\limits_{k=1}^{n+1}(x-k)}{(x-i)(-1)^{n-i+1}(i-1)!(n+1-i)!} \]

可以通过预处理的方法达到 \(\mathcal{O}(n)\) 的复杂度
CF622F The Sum of the k-th Powers
经典题目
首先我们知道自然数的 \(k\) 次幂之和一定是一个 \(k+1\) 次的多项式
可以处理出前 \(k+2\) 项之和然后直接插出来

点击查看代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<string>
#include<iostream>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099511627776,mod=1e9+7;
int n,k;
int y;
int ans;
int fl[WR],fr[WR],f[WR];
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<3)+(s<<1)+ch-'0';
        ch=getchar();
    }
    return s*w;
}
int quick_pow(int a,int b){
    int base=a,res=1;
    while(b){
        if(b&1) res=res*base%mod;
        base=base*base%mod;
        b>>=1;
    }
    return res;
}
signed main(){
    n=read(),k=read();
    f[0]=fl[0]=fr[k+3]=1;
    for(int i=1;i<=k+2;i++) fl[i]=fl[i-1]*(n-i)%mod;
    for(int i=k+2;i;i--) fr[i]=fr[i+1]*(n-i)%mod;
    for(int i=1;i<=k+2;i++) f[i]=f[i-1]*i%mod;
    for(int i=1;i<=k+2;i++){
        y=(y+quick_pow(i,k))%mod;
        int up=fl[i-1]*fr[i+1]%mod;
        int dn=f[i-1]*f[k+2-i]%mod*((k-i)&1?-1ll:1ll);
        ans=(ans+y*up%mod*quick_pow(dn,mod-2)%mod)%mod;
    }
    while(ans<0) ans+=mod;
    printf("%lld\n",ans);
    return 0;
}

P4593 [TJOI2018]教科书般的亵渎
注意到 \(n\) 的奇特的范围
显然只能对 \(m\) 搞事情，容易发现的是必然可以用 \(m+1\) 次亵渎来清场
注意到假如只有 \(1\) 个隔断点，显然只需要 \(2\) 次亵渎
贡献就是 \(\sum\limits_{i=1}^{n}i^k\) 减去隔断点的贡献
也就是说，答案可以用总体贡献减去多算的贡献
考虑对 \(a\) 排序，答案就是

\[\sum\limits_{i=0}^{m}\left(\sum\limits_{j=1}^{n-a_i}j^{m+1}-\sum\limits_{j=i+1}^{m}(a_j-a_i)^{m+1}\right) \]

然后就硬插可过

点击查看代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<string>
#include<iostream>
#include<algorithm>
#define WR WinterRain
#define int long long
using namespace std;
const int WR=1001000,INF=1099511627776,mod=1e9+7;
int m,a[WR];
int f[WR],fl[WR],fr[WR];
int read(){
    int s=0,w=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){
        if(ch=='-') w=-1;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        s=(s<<3)+(s<<1)+ch-'0';
        ch=getchar();
    }
    return s*w;
}
int quick_pow(int a,int b){
    int base=a,res=1;
    while(b){
        if(b&1) res=res*base%mod;
        base=base*base%mod;
        b>>=1;
    }
    return res;
}
int Lagrange(int x){
    fl[0]=fr[m+4]=f[0]=1;
    for(int i=1;i<=m+3;i++) fl[i]=fl[i-1]*(x-i)%mod;
    for(int i=m+3;i;i--) fr[i]=fr[i+1]*(x-i)%mod;
    for(int i=1;i<=m+3;i++) f[i]=f[i-1]*i%mod;
    int y=0,ans=0;
    for(int i=1;i<=m+3;i++){
        y=(y+quick_pow(i,m+1))%mod;
        int up=fl[i-1]*fr[i+1]%mod;
        int dn=f[i-1]*f[m+3-i]%mod*((m+1-i)&1?-1ll:1ll);
        ans=(ans+y*up%mod*quick_pow(dn,mod-2)%mod)%mod;
    }
    while(ans<0) ans+=mod;
    return ans;
}
signed main(){
    int T=read();
    while(T--){
        int n=read();m=read();
        for(int i=1;i<=m;i++) a[i]=read();
        sort(a+1,a+1+m);
        while(a[m]==n) m--,n--;
        int ans=Lagrange(n);
        for(int i=1;i<=m;i++) ans=(ans-quick_pow(a[i],m+1)+mod)%mod;
        for(int i=1;i<=m;i++) ans=(ans+Lagrange(n-a[i]))%mod;
        for(int i=1;i<=m;i++){
            for(int j=i-1;j;j--){
                ans=(ans-quick_pow(a[i]-a[j],m+1)+mod)%mod;
            }
        }
        while(ans<0) ans+=mod;
        printf("%lld\n",ans);
    }
    return 0;
}