拉格朗日插值artalter级服务

1.介绍（可忽略）

在数值分析中，拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。

2.插值

插值的定义是，

在离散数据的基础上补插连续的函数，使得这条连续函数经过所有离散数据点，这个过程就叫插值。
                                ————《百度》

通俗的讲，插值就是给你一些点，让你根据这些点确定一个多项式。

3.拉格朗日插值

定理：
    n+1个点可以确认一个n次多项式
                                    (114.514)

拉格朗日插值的原理，就是根据这一定理和手上的k+1个点，硬凑出一个k次多项式。

如果你手上有n个点，拉格朗日插值可以表示成：

/[
f(X)=/sum_{1 /le i /le n}{y_i*g(i,X)}//

g(i,X)=/prod_{1 /le j /le n,i /ne j}{/frac{X-x_j}{x_i-x_j}}
/]

这个式子看上去非常牛逼，但原理其实非常暴力:

根据上面的式子

/[g(i,X)=
/begin{cases}
1&X=i//
0&X /ne i,X/in x//
是什么都没关系的数 &X/notin x
/end{cases}
//
f(X)=
/begin{cases}
y[X]&X/in x//
用拉插算出来的f(X)&X/notin x
/end{cases}

/]

代码如下：

int fsp(int a,int b=mod-2){
	int s=1;
	while(b){
		if(b&1)s=s*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return s;
}
int f1(int n,int X,int i){
	int s1=1,s2=1;
	for(int j=1;j<=n;j++)
	{
		if(i!=j){
			s1*=(X-x[j]);
			s1%=mod;
			s2*=(x[i]-x[j]);
			s2%=mod;
		}
	}
	s2=fsp(s2);
	return s1*s2%mod;
}
int f2(int n,int X){
	int ans=0;
	for(int i=1;i<=n;i++){
		ans+=y[i]*f1(n,X,i)%mod;
		ans%=mod;
	}
	return ans;
}

拉格朗日插值复杂度是/(O(k^2)/)的

4.自然数幂和

自然数幂和，即

/[f(n)=/sum_{1/le i /le n}{i^k}
/]

是拉格朗日插值的一个重要应用。

我们知道

/[k=1时,f(n)=1+2+3+/cdots+n=/frac{n(n+1)}{2}//
k=2时，f(n)=1^2+2^2+3^2+/cdots+n^2=/frac{n(n+1)(2n+1)}{6}//
k=3时,f(n)=1^3+2^3+3^3+/cdots+n^3=(/frac{n(n+1)}{2})^2//
/]

那么，对于更大的/(k,f(n)/)有没有通用的封闭形式呢？

事实上，自然数幂和大约是没有通用的封闭形式。

但是，我们由k=1,2,31的情况可以大胆的猜一个结论：

/[
f(n)是一个k+1次多项式

/]

这个结论是对的，可以用归纳法由多项式差分证得,这里不再赘述。

那么，既然知道了这个结论，我们就可以通过算出较小情况的/(f(x)/)获得点值进行插值,/(O(k^2)/)快速(一般/(n/ggg k/))的求出/(f(n)/)了。

代码如下：

#include <bits/stdc++.h>
using namespace std;
#define int long long
const int mod = 1e9+7;
const int maxn = 2005;
int x[maxn],y[maxn];
int n;
int fsp(int x,int k = mod-2)
{
    int s=1;
    while(k)
    {
        if(k&1)s=s*x%mod;
        x=x*x%mod;
        k>>=1;
    }
    return s;
}
int fi(int n,int i,int X)
{
    int s1=1,s2=1;
    for(int j=1;j<=n;j++)
    {
        if(i!=j)
        {
            s1=s1*(X-x[j])%mod;
            s2=s2*(x[i]-x[j])%mod;
        }
    }
    return s1*fsp(s2)%mod;
}
int F(int n,int x)
{
    int s=0;
    for(int i=1;i<=n;i++)
    {
        s=s+y[i]*fi(n,i,x);
        s%=mod;
    }
    return s;
}
signed main()
{
    int n,k;
    cin>>n>>k;
    int sum=0;
    for(int i=1;i<=k+2;i++)
    {
        sum+=fsp(i,k);
        x[i]=i;
        y[i]=sum;
    }
    cout<<F(k+2,n);   
    return 0;
}