拉格朗日插值

通用模板
拉格朗日插值

给出 n + 1 n+1 n+1 x x x互不相同点 ( x i , y i ) (x_i,y_i) (xi,yi), 可以唯一确定一个 n n n次多项式

即: f ( x ) = ∑ i = 0 n a i ⋅ x i f(x)=\sum\limits_{i=0}^{n}a_i\cdotp x^i f(x)=i=0naixi 的系数 a 0 , a 1 , a 2 , ⋯   , a n a_0,a_1,a_2,\cdots,a_n a0,a1,a2,,an唯一确定

a 0 , a 1 , a 2 , ⋯   , a n a_0,a_1,a_2,\cdots,a_n a0,a1,a2,,an当作未知数, 带入 n + 1 n+1 n+1 ( x i , y i ) (x_i,y_i) (xi,yi)得到

{ a 0 ⋅ ( x 0 ) 0 + a 1 ⋅ ( x 0 ) 1 + ⋯ + a n ⋅ ( x 0 ) n = y 0 a 0 ⋅ ( x 1 ) 0 + a 1 ⋅ ( x 1 ) 1 + ⋯ + a n ⋅ ( x 1 ) n = y 1 a 0 ⋅ ( x 2 ) 0 + a 1 ⋅ ( x 2 ) 1 + ⋯ + a n ⋅ ( x 2 ) n = y 2 ⋮ a 0 ⋅ ( x n ) 0 + a 1 ⋅ ( x n ) 1 + ⋯ + a n ⋅ ( x n ) n = y n \begin{cases} a_0\cdotp (x_0)^0+a_1\cdotp (x_0)^1+\cdots+a_n\cdotp (x_0)^n=y_0\\ a_0\cdotp (x_1)^0+a_1\cdotp (x_1)^1+\cdots+a_n\cdotp (x_1)^n=y_1\\ a_0\cdotp (x_2)^0+a_1\cdotp (x_2)^1+\cdots+a_n\cdotp (x_2)^n=y_2\\ \vdots\\ a_0\cdotp(x_n)^0+a_1\cdotp (x_n)^1+\cdots+a_n\cdotp (x_n)^n=y_n\\ \end{cases} a0(x0)0+a1(x0)1++an(x0)n=y0a0(x1)0+a1(x1)1++an(x1)n=y1a0(x2)0+a1(x2)1++an(x2)n=y2a0(xn)0+a1(xn)1++an(xn)n=yn

矩阵形式:

( 1 x 0 ( x 0 ) 2 ⋯ ( x 0 ) n 1 x 1 ( x 1 ) 2 ⋯ ( x 1 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 x n ( x n ) 2 ⋯ ( x n ) n ) ( a 0 a 1 a 2 ⋮ a n ) = ( y 0 y 1 y 2 ⋮ y n ) \begin{pmatrix} 1&x_0&(x_0)^2&\cdots&(x_0)^n\\ 1&x_1&(x_1)^2&\cdots&(x_1)^n\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&x_n&(x_n)^2&\cdots&(x_n)^n \end{pmatrix} \begin{pmatrix} a_0\\a_1\\a_2\\ \vdots\\a_n \end{pmatrix}= \begin{pmatrix} y_0\\y_1\\y_2\\ \vdots\\y_n \end{pmatrix} 111x0x1xn(x0)2(x1)2(xn)2(x0)n(x1)n(xn)na0a1a2an=y0y1y2yn

用高斯消元可以在 O ( n 3 ) O(n^3) O(n3)内得到 a 0 , a 1 , a 2 , ⋯   , a n a_0,a_1,a_2,\cdots,a_n a0,a1,a2,,an

构造多项式: l i ( x ) = { 1 ( x = x i ) 0 ( x = x j ⋂ i ≠ j ) l_i(x)=\begin{cases}1\qquad(x=x_i)\\ 0\qquad(x=x_j \bigcap i\neq j) \end{cases} li(x)={1(x=xi)0(x=xji̸=j)

**那么: ** f ( x ) = ∑ i = 0 n l i ( x ) ⋅ y i f(x)=\sum\limits_{i=0}^nl_i(x)\cdotp y_i f(x)=i=0nli(x)yi

令: l i ( x ) = ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x i − 1 ) ( x − x i + 1 ) ⋯ ( x − x n ) ( x i − x 0 ) ( x i − x 1 ) ⋯ ( x i − x i − 1 ) ( x i − x i + 1 ) ⋯ ( x i − x n ) = ∏ j = 0 ⋂ j ≠ i n ( x − x j ) ∏ j = 0 ⋂ j ≠ i n ( x i − x j ) l_i(x)=\frac{(x-x_0)(x-x_1)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)(x_i-x_1)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}=\frac{\prod\limits_{j=0\bigcap j\neq i}^n (x-x_j)}{\prod\limits_{j=0\bigcap j\neq i}^n(x_i-x_j)} li(x)=(xix0)(xix1)(xixi1)(xixi+1)(xixn)(xx0)(xx1)(xxi1)(xxi+1)(xxn)=j=0j̸=in(xixj)j=0j̸=in(xxj)

//(x[0],y[0]),……,(x[n],y[n])    n+1个点
//插出n次多项式, 返回带入x0的值
ll lagrange(int x[],int y[],int n,int x0)
{
    ll res=0;
    for(int i=0;i<=n;++i)
    {
        ll fz=1,fm=1;
        for(int j=0;j<=n;++j)
        {
            if(j==i)continue;
            fz=md(fz*(x0-x[j]));
            fm=md(fm*(x[i]-x[j]));
        }
        (res+=y[i]*fz%mod*qpow(fm,mod-2)%mod)%=mod;
    }
    return res;
}
自然幂和

求: ∑ i = 1 n i k \sum\limits_{i=1}^ni^k i=1nik

//快速幂

const int __=1e5+5;
int pk[__+5],num[__+5],sum[__+5];

//线性筛: pk[i]=qpow(i,k)%mod
void powk(int k)
{
    num[0]=0,pk[1]=sum[1]=1;
    for(int i=2;i<=k+1;++i)pk[i]=0;
    for(int i=2;i<=k+1;++i)
    {
        if(!pk[i])
        {
            num[++num[0]]=i;
            pk[i]=qpow(i,k);
        }
        for(int j=1;j<=num[0] && i*num[j]<=__;++j)
        {
            int &p=num[j];
            pk[i*p]=1ll*pk[i]*pk[p]%mod;
            if(!(i%p))break;
        }
        sum[i]=(sum[i-1]+pk[i])%mod;
    }
}

//预处理: fac[]阶乘, inv[]逆元, facinv[]阶乘逆元

//(x[i]=x+i,y[i])    n+1个x连续的点
//插出n次多项式, 返回带入x0的值
ll pre[__],saf[__];

ll lagrange(int x,int y[],int n,int x0)
{
    pre[0]=x0-x,saf[n]=x0-x-n;
    for(int i=1;i<=n;++i)
    {
        pre[i]=md(pre[i-1]*(x0-x-i));
        saf[n-i]=md(saf[n-i+1]*(x0-x-n+i));
    }
    ll res=0;
    for(int i=0;i<=n;++i)
    {
        ll fz=1;
        if(i!=0)fz=md(fz*pre[i-1]);
        if(i!=n)fz=md(fz*saf[i+1]);
        ll fm=md(facinv[i]*facinv[n-i]);
        if((n-i)&1)fm=mod-fm;
        (res+=md(y[i]*md(fz*fm)))%=mod;
    }
    return res;
}

int main()
{
    init();
    int _;for(sf("%d",&_);_;--_)
    {
        ll n;int k;
        sf("%lld%d",&n,&k);
        powk(k);
        pf("%lld\n",lagrange(0,sum,k+1,n%mod));
    }
    return 0;
}

复杂度: O ( q ⋅ k ) O(q\cdot k) O(qk)

  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值