拉格朗日插值与牛顿插值

2016.03.23 18:45 Wed | 150次阅读 | 旧日oi | 固定链接 | 源码

前言

今天上午又考了一道插值的题,结果n^3的高斯消元就能过,感觉不太优美,下午改完题之后又去看了看插值的东西,为了防止再忘掉,写篇文章总结一下。
我并不会什么黑科技,所以只总结下两种n^2的插值法。

插值是什么

给出平面上一些点,要求搞出一个多项式函数来拟合这些点,可能问这个多项式是什么,也可能问这个多项式在某个点的值是什么,本质上是相同的。
常见的插值方法有两种:牛顿插值和拉格朗日插值。

牛顿插值

网上有篇文章是从公式出发,然后引出差商的,并且也并没有介绍差商。我觉得不太好理解,差商真是一个很奇妙的东西。
定义f[x[i],x[j]]表示x[i]到x[j]的差商,x[i],f[x[i]]表示x=x[i]时的函数值是f[x[i]]。

假设现在要求拟合出的多项式f(x)。

然后这么一直迭代搞下去就行了,差商f[x[0],x[n+1]]之后就都是0了,那么这个公式最多也就n次。
依次带回原式就可以得到一个很优美的公式

感觉差商和函数值长的差不多?注意区分下就好了……
我们可以n^2预处理出差商表,然后如果要求多项式f[x]在某个点的值的时候直接带入就好了。
如果要求整个多项式就要每次多项式乘法了,好在乘法是可以O(n)完成的,所以求一个多项式的复杂度也是n^2.

拉格朗日插值

这个我先给公式:

怎么验证呢,可以把x0...xn带入就行了,带入x_i时,除了li(x)=1之外,其它都是0……
求某个点的值的话,只要维护一个前缀积和后缀积就好了。
求多项式的值的话
我们可以先把$(x-x_0)...(x-x_n)$求出来,每次再模一个一次多项式就可以了,复杂度也是n^2。

题目

hdu5275

题目大意

给出n个点,Q个询问,问用[l,r]之间的点拟合出的多项式在某个点的值。

题解

先考虑牛顿插值,直接把f[x[l],x[l+1]]....f[x[l],x[r]]带入就好了,相当简单。
然后是拉格朗日插值。
可以不考虑不乘的那项,然后预处理分母,到时候乘个逆元就好了。
如果x恰好和某个给出的点相同,直接判一下就好了。
否则也是都乘起来每次求逆元。
两种办法都是预处理n^2,每次询问O(n)

mycode

牛顿插值

#include <bits/stdc++.h>
using namespace std;
#define N 250005
#define ll long long
#define mod 1000000007
int n,Q,l,r,q;
ll x[N],y[N],inv[N];
ll f[3005][3005];
ll get_inv(int x){return x<0?-inv[-x]:inv[x];}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("tt.in","r",stdin);
    #endif
    cin>>n;
    for(int i=1;i<=n;i++) scanf("%lld%lld",x+i,y+i);
    inv[1]=1;for(int i=2;i<=250000;i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;    
    for(int i=1;i<=n;i++) f[0][i]=y[i];
    for(int i=1;i<=n;i++)
        for(int j=1;j+i<=n;j++)
            f[i][j]=((f[i-1][j+1]-f[i-1][j])*get_inv(x[i+j]-x[j])%mod+mod)%mod;
    for(cin>>Q;Q--;)
    {
        scanf("%d%d%d",&l,&r,&q);
        ll ans=0,now=1;
        for(int i=0;i<=r-l;i++)
        {
            ans=(ans+f[i][l]*now)%mod;
            now=(now*(q-x[l+i])%mod+mod)%mod;
        }
        printf("%lld\n",ans);
    }
}

拉格朗日插值

#include <bits/stdc++.h>
using namespace std;
#define N 250005
#define ll long long
#define mod 1000000007
int n,Q,l,r,q;
ll x[N],y[N],inv[N],_inv[N];
ll f[3005][3005];
ll get_inv(int x){return x<0?-inv[-x]:inv[x];}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("tt.in","r",stdin);
    #endif
    cin>>n;
    for(int i=1;i<=n;i++) scanf("%lld%lld",x+i,y+i);
    inv[1]=1;
    for(int i=2;i<=250000;i++) 
    {
        inv[i]=(mod-mod/i)*inv[mod%i]%mod;    
        _inv[i]=_inv[i-1]*inv[i];
    }
    for(int i=1;i<=n;i++) f[0][i]=y[i];
    for(int i=1;i<=n;i++)
    {
        f[i][i]=1;
        for(int j=i-1;j;j--) f[i][j]=f[i][j+1]*get_inv(x[i]-x[j])%mod;
        for(int j=i+1;j<=n;j++) f[i][j]=f[i][j-1]*get_inv(x[i]-x[j])%mod;
    }
    for(cin>>Q;Q--;)
    {
        scanf("%d%d%d",&l,&r,&q);
        ll ans=0,tmp=1;int pos=0;
        for(int i=l;i<=r;i++)
        {
            if(x[i]!=q) tmp=tmp*(q-x[i])%mod;
            else pos=i;
        }
        if(pos) ans=f[pos][l]*f[pos][r]%mod*y[pos]%mod*tmp%mod;
        else for(int i=l;i<=r;i++) 
            ans=(ans+y[i]*f[i][l]%mod*f[i][r]%mod*tmp%mod*get_inv(q-x[i])%mod)%mod;
        printf("%lld\n",(ans+mod)%mod);
    }
}

CF622F

题目大意

求1到n的k次方和,n<=1e9,k<=1e6

解题思路

1....n的k次方和是关于n的k次多项式,那么我们可以只用k+1个点把这个多项式拟合出来再把n带入就好了。
而朴素的插值是n^2的,我们可以选一些特殊点1,2....k+1,用拉格朗日插值法做到nlogn
问题主要在维护分母部分,注意到从x[i]=i转移到x[i+1]=i+1并不会有很多东西变化,只有一次乘法和一次除法,推推公式,费马小定理求逆元,就可以轻松的O(logn)转移了。

mycode

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define mod 1000000007
int n,k;
ll ans[1000005];
ll qpow(ll x,ll k)
{
    ll ret=1;
    while(k)
    {
        if(k&1) ret=ret*x%mod;
        k>>=1;x=x*x%mod;
    }
    return ret;
}
int main()
{
    #ifndef ONLINE_JUDGE
    freopen("tt.in","r",stdin);
    #endif
    cin>>n>>k;
    for(int i=1;i<=k+2;i++) ans[i]=(ans[i-1]+qpow(i,k))%mod;
    if(n<=k+2) 
    {
        printf("%I64d\n",ans[n]);
        return 0;
    }
    ll tmp=1,fin=0,d=1;
    for(int i=2;i<=k+2;i++) d=d*(1-i+mod)%mod; 
    for(int i=1;i<=k+2;i++) tmp=tmp*(n-i)%mod;
    for(int i=1;i<=k+2;i++) 
    {
        fin=(fin+ans[i]*tmp%mod*qpow(n-i,mod-2)%mod*qpow(d,mod-2))%mod;
        d=d*qpow(i-k-2+mod,mod-2)%mod*i%mod;
    }
    printf("%I64d\n",fin);
}

TCO2014 Round 3B Div1 1000

题目大意

给出一个n个点的树,允许修至多改k条边的连法,使得修改完还是一棵树,问能得到多少种不同的树。

题解

我们可以把原有的边权值设为1,非原有的边边权设为x,通过矩阵树定理可以得到所有可能的生成树之积,注意到这个积是一个由x构成的n-1次多项式,那么我们令x=0...n进行矩阵树定理,就可以得到n+1个对应的值。
有了n个点之后,我们就可以利用牛顿插值轻松的求出这个多项式了。

mycode

#include <bits/stdc++.h>
using namespace std;
#define mod 998244353
#define ll long long
int n,k,cnt,f[55];
int pos[55][55];
ll a[55][55],c[55],d[55],x[55],y[55],ans[55];
ll b[55][55];
ll qpow(ll x,ll k)
{
    ll ret=1;
    while(k)
    {
        if(k&1) ret=ret*x%mod;
        k>>=1;x=x*x%mod;
    }
    return ret;
}
ll gauss(int equ,int var)
{
    int i,j,k;
    ll ret=1;
    for(i=0;i<equ;i++)
    {
        if(!a[i][i]) return 0;
        ret=ret*a[i][i]%mod;
        for(j=i;j<equ;j++)
        if(j!=i&&a[j][i])
        {
            ll tmp=qpow(a[i][i],mod-2)*a[j][i]%mod;
            for(k=i;k<=var;k++) 
                a[j][k]=(a[j][k]-a[i][k]*tmp%mod+mod)%mod;
        }
    }
    return ret;
}
int main()
{
    freopen("tree.in","r",stdin);
    freopen("tree.out","w",stdout);
    cin>>n>>k;
    for(int i=1;i<n;i++) 
    {
        scanf("%d",&f[i]);
        pos[f[i]][i]=pos[i][f[i]]=1;
    }
    for(int val=0;val<n;val++)
    {
        memset(a,0,sizeof(a));
        for(int i=0;i<n;i++) 
        for(int j=0;j<n;j++) 
        if(i!=j)
        {
            if(pos[i][j]) a[i][j]=mod-1,a[i][i]++;
            else a[i][j]=mod-val,a[i][i]+=val;
        }
        x[val]=val;y[val]=gauss(n-1,n-1);
    }
    for(int i=0;i<n;i++) b[0][i]=y[i];
    for(int j=1;j<n;j++)
    {
        for(int i=0;i+j<n;i++)
        b[j][i]=(b[j-1][i]-b[j-1][i+1]+mod)%mod*qpow((x[i]-x[i+j]+mod)%mod,mod-2)%mod;
    }
    ans[0]=b[0][0];c[0]=1;
    for(int i=1;i<n;i++)
    {
        memset(d,0,sizeof(d));
        for(int j=0;j<n;j++) 
        {
            d[j+1]+=c[j];
            d[j]-=x[i-1]*c[j];
        }
        for(int j=0;j<n;j++) c[j]=(d[j]%mod+mod)%mod;
        for(int j=0;j<n;j++) (ans[j]+=b[i][0]*c[j])%=mod;
    }
    int ret=0;
    for(int i=0;i<=k;i++) ret=(ret+ans[i])%mod;
    printf("%lld\n",ret);
}

对比

拉格朗日插值的优势:通过某些特殊点的值可以线性推出任意点的值。
拉格朗日插值的劣势:代码复杂度略大,子串插值难度大。
牛顿插值的优势:好写,子串插值方便
牛顿插值的劣势:复杂度过于稳定,没什么优化空间。