拉格朗日插值法

前言

在数值分析中，拉格朗日插值法是一种多项式插值方法。

拉格朗日插值法

下文，本蒟蒻可能有说的不妥甚至不对的地方， 欢迎大神来打脸。
定义多项式函数$f(x)$，则点值运算就是给定$x$，求$y=f(x)$
插值是点值的逆运算，就是给定$n$次多项式的$n+1$个点值表达$(x_i,y_i)$，要求出$f(x)$各次项系数
插值法有许多种，当然有名的有拉格朗日插值法与牛顿插值法。
有一个结论，n次的函数可以由n+1个点确定。例如两个点确定一次函数；三个点确定二次函数。因为待定系数包括常数项有n+1项，所以要有n+1个方程来解。

拉格朗日插值基函数

对于每个$(x_k,y_k)$，定义

$$l_k(x)=\prod_{i=0,i\neq k}^n\dfrac{x-x_i}{x_k-x_i}$$
这个函数有一个好处，就是 $\forall_{i\leq n}$，当且仅当 $i=k$时， $l_k(x_i)=1$，否则 $l_k(x_i)=0$
于是把每个点值处的基函数求和
$$L(x)=\sum_{i=0}^n y_il_i(x)$$
就能符合每个点值了， $L(x)$即为所求多项式函数

应用

自然数幂和

引理：$f^k(n)=$$\sum_{i=1}^ni^k$必定能用$k+1$次多项式函数表示

我们只要取$k+2$个点值表达，即可插值出这个多项式函数
为了简便直接取$(0,f^k(0)),(1,f^k(1))...(k+1,f^k(k+1))$
代入$L(x)$，得

$$f^k(n)=L(n)=\sum_{j=0}^{k+1}\frac{(-1)^{k+1-j}f^k(j)n(n-1)...(n-j+1)(n-j-1)...(n-k-1)}{j!(k+1-j)!}$$
但是它的局限性在于，需要预处理出阶乘的逆元，因此对模数有要求
时间复杂度可以做到 $O(k)$

例题 51nod1258 Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define fd(i,b,a) for(ll i=b;i>=a;i--)
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define mset(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
const ll mo=1e9+7,K=5e4+10;
ll k,fk[K],fac[K],inv[K],pre[K],suf[K],pr[K];
bool bz[K];
ll qmi(ll x,ll n)
{
    ll t=1;
    for(x%=mo;n;n>>=1,x=x*x%mo)
        if(n&1) t=t*x%mo;
    return t;
}
void sieve(ll n)
{
    mset(fk,0);mset(inv,0);mset(fac,0);
    fk[0]=0,fk[1]=fac[0]=inv[0]=fac[1]=inv[1]=1;
    pr[0]=0;
    fo(i,2,n)
    {
        fac[i]=fac[i-1]*i%mo,fk[i]=(fk[i]+fk[i-1])%mo;
        if(!bz[i]) pr[++pr[0]]=i,fk[i]=(fk[i-1]+qmi(i,k))%mo,inv[i]=qmi(i,mo-2);
        inv[i]=inv[i]*inv[i-1]%mo;
        fo(j,1,pr[0])
        {
            ll x=i*pr[j];
            if(x>n) break;
            bz[x]=1;
            inv[x]=fac[i-1]*inv[i]%mo*fac[pr[j]-1]%mo*inv[pr[j]]%mo;
            fk[x]=(fk[x]+(fk[i]-fk[i-1]+mo)*(fk[pr[j]]-fk[pr[j]-1]+mo)%mo)%mo;
            if(i%pr[j]==0) break;
        }
    }
}
ll calc(ll n)
{
    if(!n) return 0;
    if(!k) return n;
    if(n<=k+1) return fk[n];
    ll ans=0;
    int t=0;
    pre[0]=1;
    fo(i,n-k-1,n) ++t,pre[t]=pre[t-1]*(i%mo)%mo;
    t=0;
    suf[0]=1;
    fd(i,n,n-k-1) ++t,suf[t]=suf[t-1]*(i%mo)%mo;
    fo(j,0,k+1)
    {
        ll t=fk[j]*pre[k-j+1]%mo*suf[j]%mo*inv[j]%mo*inv[k+1-j]%mo;
        if((k+1-j)&1) ans=(ans-t+mo)%mo;
        else ans=(ans+t)%mo;
    }
    return ans;
}
int main()
{
    int T;ll n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld %lld",&n,&k);
        sieve(k+1);
        printf("%lld\n",calc(n));
    }
    return 0;
}