拉格朗日插值artalter级服务

拉格朗日插值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)\)的

P4781 【模板】拉格朗日插值

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;
}

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{aligned} &\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 \dots \times (i-1))\times((n-i)\times(n-i-1)\times(n-i-2)\dots \times 1)}\\ =&(-1)^{n-i}\frac{1}{(i-1)!(n-i)!}\\ \end{aligned} \]

于是，我们预处理出分母，就可以在\(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

结束

武藏老婆镇楼（诶嘿终于贴了一张我抽到的老婆了）