拉格朗日插值artalter级服务

拉格朗日插值artalter级服务

1.介绍（可忽略）

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

2.插值

插值的定义是，

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

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

3.拉格朗日插值

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

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

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

$$\begin{gathered} f(X)=\sum _{1\le i\le n}{y}_i\ast g(i,X) \\ g(i,X)=\prod _{1\le j\le n,i\ne j}\frac{X-x_j}{x_i-x_j} \end{gathered}$$
这个式子看上去非常牛逼，但原理其实非常暴力:

根据上面的式子

$$\begin{gathered} g(i,X)=\{\begin{array}{ll}1& X=i\\ 0& X\ne i,X\in x\\ 是什么都没关系的数& X\notin x\end{array} \\ f(X)=\{\begin{array}{ll}y[X]& X\in x\\ 用拉插算出来的f(X)& X\notin x\end{array} \end{gathered}$$
代码如下：

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)$的

P4781 【模板】拉格朗日插值

4.自然数幂和

自然数幂和，即
$$f(n)=\sum _{1\le i\le n}{i}^k$$
是拉格朗日插值的一个重要应用。

我们知道
$$\begin{gathered} 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 \end{gathered}$$

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

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

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

$$

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

$$

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

那么，既然知道了这个结论，我们就可以通过算出较小情况的$f(x)$获得点值进行插值,$O(k^2)$快速(一般$n⋙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;
}

5.优化

我们发现，如果拉格朗日插值中的点值可以自选，我们就可以对它进行优化。

如果我们选的$x$是从$1$开始连续的,即$x_i=i,y_i=f(i)$,我们可以预处理出
$$
MUL=\prod_{1\le i \le n}(X-x_i)

$$
$g(X,i)$的分子变成了$MUL\times (X-x_i{)}^{-1}$

而分母同样可以优化：
$$\begin{array}{rl}& \frac{1}{{\prod }_{1\le j\le n,j\ne i}(x_i-x_j)}\\ =& \frac{1}{{\prod }_{1\le j<i}(i-j){\prod }_{i<j\le n}(i-j)}\\ =& (-1{)}^{n-i}\frac{1}{{\prod }_{1\le j<i}(i-j){\prod }_{i<j\le n}(j-i)}\\ =& (-1{)}^{n-i}\frac{1}{(1\times 2\times 3\times \cdots \times (i-1))\times ((n-i)\times (n-i-1)\times (n-i-2)\cdots \times 1)}\\ =& (-1{)}^{n-i}\frac{1}{(i-1)!(n-i)!}\end{array}$$

于是，我们预处理出分母，就可以在$O(\log MOD)$的复杂度内求出g。整个拉插的复杂度也就为$O(n\log MOD)$

代码如下：

int x[maxn],y[maxn],fr[maxn],FR,frac[maxn],fa[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=FR*fsp(X-x[i])%mod;
    return s1*fr[i]%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;
}
void init()
{
    frac[0]=1;
    for(int i=1;i<=maxn-5;i++)frac[i]=frac[i-1]*i%mod;
}
int numpow(int n,int k)
{
    memset(x,0,sizeof(x));
    memset(y,0,sizeof(y));
    memset(fr,0,sizeof(fr));
    int sum=0;
    for(int i=1;i<=k+2;i++)
    {
        sum+=fsp(i,k);
        sum%=mod;
        x[i]=i;
        y[i]=sum;
    }
    FR=1;
    for(int i=1;i<=k+2;i++)
    {
        fr[i]=((k - i) & 1 ? -1 : 1)*fsp(frac[i-1])%mod*fsp(frac[k+2-i])%mod;
        FR*=(n-i);FR%=mod;
    }
    if(n<=k+2)
        return y[n];
    return (F(k+2,n)+mod)%mod; 
}

板子题地址：CF622F The Sum of the k-th Powers