浅谈拉格朗日插值
在数值分析中,拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律,而不少函数都只能通过实验和观测来了解。拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日(插值)多项式。
——百度百科
通俗地说,拉格朗日插值法可以找出一个恰好经过直角坐标系内(n)个给定点的函数
众所周知,(n)个点能唯一确定的多项式最高次数是(n-1)次
这个可以两点确定一次函数,三点确定二次函数来推出,或者我们由方程组有唯一解的充要条件也能得出
知道了这个之后我们就容易想到直接用高斯消元来搞出多项式的系数,但是复杂度消耗太大,是(O(n^3))
而拉格朗日插值法就是一种一般可以在(O(n^2))的复杂度下求出多项式的方法(不过通常用来求一个点值,以下的讲述一般以此为参考)
一般的拉格朗日插值法
我们设要求的多项式为(f(x)),点的坐标为((x_i,y_i)),我们要找多项式在(x_0)处的取值
先上公式,由拉格朗日插值法:
看起来不是很好理解,其实很简单,我们把原来的某个给定点(x_k)代入以下有:
容易发现,当(k ot=i)时,后面的(frac{x_k-x_j}{x_i-x_j})的分子总有一项是(0),此时(prod_{i ot= j} frac{x_k-x_j}{x_i-x_j}=0)
当(k=i)时,后面的(frac{x_k-x_j}{x_i-x_j})上下完全相同,此时(prod_{i ot= j} frac{x_k-x_j}{x_i-x_j}=1)
即对于(f(x_k))来说,这个多项式的确给出了对应的(y_k)的值
不难发现这个方法对所有点都适用,因此它是正确的
从上面的式子可以看出每次计算要枚举两次,因此复杂度很简单,就是(O(n^2))
在(x)取值连续时的插值法
因为很多时候我们做题都是先发现某个函数是多少次的多项式,然后自己随意取一些值代入插值
这样的话为了省事横坐标的取值完全可以从(1)开始连续取,那么我们把上面的式子中的(x_i)换成(i)就有:
考虑怎么快速求(prod_{i ot= j} frac{k-j}{i-j}),我们分别考虑:
分子的话容易发现就是(frac{prod_{t=1}^n x_0-t}{x_0-i})
分母比较复杂,(i-j)的累乘可以分成两个阶乘部分,因此推导一下就是((-1)^{n-i}cdot i!cdot(n-i)!)
这样我们一般就可以(O(nlog n))来算了((log)主要是有求逆元的过程,当然你预处理(O(n))求也没有问题)
重心拉格朗日插值法
我们考虑到朴素的拉格朗日插值每次多加入一个点时就要整个重新算过,很浪费时间
那么能不能把重复算的一些东西利用起来?
我们对于
把分子提取出来,设为(g=prod_{i=1}^n x_0-x_i),则此时:
设一个(t_i=frac{y_i}{prod_{i ot= j} x_i-x_j}),则:
因此每次多加入一个点只需要重新(O(n))算它的(t_i)就好了
拉格朗日插值的应用以及常用解题思路
一个经典例子:自然数幂和,即求(sum_{i=1}^n i^k)
之前也提到过用第二类斯特林数做的方法,但那种方法是(O(k^2))的,不够优秀
但是现在我们观察这个式子,如果不看求和的话(i^k)就是一个(k)次多项式
那么前缀和(其实就是差分)之后,次数要(+1),即此时答案为一个关于(n)的(k+1)多项式
那我们直接代入(k)个值(取值连续,反正自己定)之后插值算即可,复杂度(O(klog k))
那么很多具体题目怎么办呢,一般就是先推出某个式子,然后证明它是关于某个自变量的多少次函数(一定要判断出次数!),然后自己选一些点(一般是连续的)代入插值即可
其实做了一些题目之后就很套路了
拉格朗日插值的模板
模板看Luogu P4781 【模板】拉格朗日插值 ,就是用一般的(O(n^2))来做就好了
CODE
#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=2005,mod=998244353;
int n,x[N],y[N],k;
inline int sub(CI x,CI y)
{
int t=x-y; return t<0?t+mod:t;
}
inline int inv(int x,int p=mod-2,int mul=1)
{
for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline int Lagrange(CI n,int *x,int *y,CI k)
{
int ret=0; for (RI i=1;i<=n;++i)
{
int s1=1,s2=1; for (RI j=1;j<=n;++j) if (i!=j)
s1=1LL*s1*sub(k,x[j])%mod,s2=1LL*s2*sub(x[i],x[j])%mod;
(ret+=1LL*y[i]*s1%mod*inv(s2)%mod)%=mod;
}
return ret;
}
int main()
{
scanf("%d%d",&n,&k); for (RI i=1;i<=n;++i) scanf("%d%d",&x[i],&y[i]);
return printf("%d",Lagrange(n,x,y,k)),0;
}
一些入门例题
- CF622F The Sum of the k-th Powers 提到过的自然数幂和
- Luogu P3270 [JLOI2016]成绩比较 需要容斥+组合数学+广义容斥化式子的自然数幂和
- Luogu P4463 [国家集训队] calc DP+多项式差分系数推导之后用拉格朗日插值算答案