【拉格朗日插值法】

最新推荐文章于 2022-09-24 19:57:24 发布

原创最新推荐文章于 2022-09-24 19:57:24 发布 · 313 阅读

0 ·

CC 4.0 BY-SA版权

参考博客：https://www.cnblogs.com/ECJTUACM-873284962/p/6833391.html
定义

对某个多项式函数，已知有给定的k + 1个取值点：

$(x_{0},y_{0}),\ldots ,(x_{k},y_{k})$

其中 $x_{j}$ 对应着自变量的位置，而 $y_{j}$ 对应着函数在这个位置的取值。

假设任意两个不同的 $x_j$ 都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：

$L(x):=\sum _{{j=0}}^{{k}}y_{j}\ell _{j}(x)$

其中每个 $\ell _{j}(x)$ 为拉格朗日基本多项式（或称插值基函数），其表达式为：

$\ell _{j}(x):=\prod _{{i=0,\,i\neq j}}^{{k}}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{{j-1}})}{(x_{j}-x_{{j-1}})}}{\frac {(x-x_{{j+1}})}{(x_{j}-x_{{j+1}})}}\cdots {\frac {(x-x_{{k}})}{(x_{j}-x_{{k}})}}$

拉格朗日基本多项式 $\ell _{j}(x)$ 的特点是在 $x_{j}$ 上取值为1，在其它的点 $x_{i},\,i\neq j$ 上取值为0。

代码

//1.先调用 init（n）初始化 逆元数组。
//2.估计出答案的最高次数，然后构造一个比最高次数多一项的插基数组。
//3.然后调用polyval求第n项。
typedef long long LL;
const int MOD=1e9+7;
const int N=1020;
LL pod(LL x,LL n)
{
    LL ret=1;
    while(n)
    {
        if(n&1)ret=ret*x%MOD;
        n>>=1;
        x=x*x%MOD;
    }
    return ret;
}
namespace Polyval
{
LL fac[N],invv[N],p1[N],p2[N];

void init(int n)
{
    fac[0]=fac[1]=invv[0]=invv[1]=1;
    for(int i=2; i<=n; i++)
        fac[i]=fac[i-1]*i%MOD;
    invv[n]=pod(fac[n],MOD-2);
    for(int i=n; i>1; i--)
        invv[i-1]=invv[i]*i%MOD;
}
LL polyval(int d,LL *a,LL n)
{
    if(n<=d)return a[n];
    p1[0]=1;
    p2[d]=1;
    for(LL i=0; i<=d; i++)
        p1[i+1]=p1[i]*(n-i)%MOD;
    for(LL i=d; i>0; i--)
        p2[i-1]=p2[i]*(n-i)%MOD;
    LL ans=0;
    for(int i=0; i<=d; i++)
    {
        LL tem=a[i]*p1[i]%MOD*p2[i]%MOD*invv[i]%MOD*invv[d-i]%MOD;
        cout<<tem<<endl;
        if((d-i)&1)ans=(ans-tem+MOD)%MOD;
        else ans=(ans+tem)%MOD;
    }
    ans=(ans+MOD)%MOD;
    return ans;
}
}
LL b[N],a[N];
//求k次方的前n项和
int main()
{
    int n,k;
    Polyval::init(1010);
    while(scanf("%d%d",&n,&k)!=EOF)
    {
        b[0]=0;
        for(LL i=1;i<=k+1;i++)  //和最高为k+1项
        {
            b[i]=(pod(i,k)+b[i-1])%MOD; //构造和的k+2个点
        }
        printf("%lld\n",Polyval::polyval(k+1,b,n));
    }
}

杜教模板