拉格朗日插值法学习笔记

STEP1 问题引入

现在给定在平面直角坐标系中的 k + 1 k+1 k+1 个点,记 P i ( x i , y i ) P_i(x_i,y_i) Pi(xi,yi)。现在要求一个 k k k 次多项式函数 y = ∑ i = 0 k a i x i y=\sum_{i=0}^ka_ix^i y=i=0kaixi,使得它经过这 k + 1 k+1 k+1 个点,即 ∀ i ∣ 0 ≤ i ≤ k , y ∣ x = x i = y i \forall_{i|0\le i\le k},y|_{x=x_i}=y_i i0ik,yx=xi=yi。输出 x = n x=n x=n y y y 的取值。

我们首先想到的是把 k + 1 k+1 k+1 个点的坐标分别代入 y = ∑ i = 0 k a i x i y=\sum_{i=0}^ka_ix^i y=i=0kaixi,则可以得到 k + 1 k+1 k+1 个方程。于是便可以用高斯消元法求出所有的 a i a_i ai,最后把 n n n 代进去算就好了时间复杂度为 O ( k 3 ) O(k^3) O(k3)。这种方法效率非常低,需要用更高效的方法代替。

STEP2 拉格朗日插值公式

存在性

先考察多项式 l i = ∏ j ≠ i x − x j x i − x j l_i=\prod_{j\not=i}\frac{x-x_j}{x_i-x_j} li=j=ixixjxxj。不难发现,如果 x = x i x=x_i x=xi,那么每个分数分子与分母相等,即每个分数都为 1 1 1,则此时 l i = 1 l_i=1 li=1;如果 x ≠ x i x\not=x_i x=xi x x x 为某个点 P j P_j Pj 的横坐标,那么分子有一项为 x − x j = 0 x-x_j=0 xxj=0,那么此时 l i = 0 l_i=0 li=0。像这样的多项式,我们称之为拉格朗日基本多项式
那么如果将 l i l_i li 乘上 y i y_i yi ,就可以得到在 x i x_i xi 上取值为 y i y_i yi,在 x j x_j xj j ≠ i j\not=i j=i)上取值为 0 0 0 的一系列多项式。我们把所有这种多项式累加起来,就可以得到要求的多项式 L ( x ) = ∑ i = 0 k y i l i = ∑ i = 0 k y i ∑ 0 ≤ j ≤ k , j ≠ i x − x j x i − x j L(x)=\sum_{i=0}^ky_il_i=\sum_{i=0}^ky_i\sum_{0\le j\le k,j\not=i}\frac{x-x_j}{x_i-x_j} L(x)=i=0kyili=i=0kyi0jk,j=ixixjxxj
时间复杂度是 O ( k 2 ) O(k^2) O(k2)
由于每个点的横坐标不相等, i ≠ j i\not=j i=j,故 x i ≠ x j x_i\not=x_j xi=xj,即 x i − x j ≠ 0 x_i-x_j\not=0 xixj=0。也就是说,上式中的所有分母不等于零,故这样的多项式是存在的。

唯一性

假设有两个不同的 k k k 次多项式 L 1 ≠ L 2 L_1\not=L_2 L1=L2,记 L = L 1 − L 2 L=L_1-L_2 L=L1L2,则 L L L 的次数不大于 k k k。又 L 1 L_1 L1 L 2 L_2 L2 x 0 , x 1 , … , x k x_0,x_1,\dots,x_k x0,x1,,xk 上的取值相同,则 L L L 在这些点上取值都为零,也就是说 L = a ∏ i = 0 k ( x − x i ) L=a\prod_{i=0}^k(x-x_i) L=ai=0k(xxi)。若 a ≠ 0 a\not=0 a=0,则这是一个 k + 1 k+1 k+1次多项式,与它次数不大于 k k k 相矛盾,故 a = 0 a=0 a=0,那么 L = 0 L=0 L=0 L 1 = L 2 L_1=L_2 L1=L2。这又与 L 1 ≠ L 2 L_1\not=L_2 L1=L2 相矛盾。故有且只有一个 k k k 次多项式经过平面上的 k + 1 k+1 k+1 个点。

代码模板

#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
int n;
ll x[2001],y[2001],mod=1e9+7;
ll k,ans,w,l;
inline ll pw(ll a,int b){  //快速幂
	if(b==0) return 1ll;
	if(b==1) return a;
	ll res=pw(a,b>>1);
	(res*=res)%=mod;
	if(b&1) (res*=a)%=mod;
	return res;
}
inline ll sub(ll a,ll b){  //模下的减法
	return (a-b>0)?(a-b):(a-b+mod);
}
inline ll inv(ll a){  //乘法逆元
	return pw(a,998244351);
}
inline void lagrange(){
	for(register int i=0;i<=k;i+=1){
		l=1ll; w=1ll;
		for(register int j=0;j<=k;j+=1){
			if(i==j) continue;
			(l*=sub(n,x[j]))%=mod;  //计算分子
			(w*=sub(x[i],x[j]))%=mod;  //计算分母
		}
		(l*=inv(w))%=mod;  //拉格朗日基本多项式的值
		(ans+=l*y[i])%=mod;  //累加答案
	}
	return;
}
int main(){
	scanf("%d%lld",&n,&k);
	for(register int i=0;i<=k;i+=1){  //下标从0开始
		scanf("%lld%lld",&x[i],&y[i]);
	}
	lagrange();
	printf("%lld\n",ans);
	return 0;
}

STEP3 算法优化

重心拉格朗日插值

定义 l ( x ) = ∏ i = 0 k ( x − x i ) l(x)=\prod_{i=0}^k(x-x_i) l(x)=i=0k(xxi),那么有 l i ( x ) = l x − x i l_i(x)=\frac{l}{x-x_i} li(x)=xxil
原式变成 L ( x ) = ∑ i = 0 k y i l ( x ) ( x − x i ) ∏ 0 ≤ j ≤ k , j ≠ i ( x i − x j ) L(x)=\sum_{i=0}^k\frac{y_il(x)}{(x-x_i)\prod_{0\le j\le k,j\not=i}(x_i-x_j)} L(x)=i=0k(xxi)0jk,j=i(xixj)yil(x)。定义重心权 w i = 1 ∏ 0 ≤ j ≤ k , j ≠ i ( x i − x j ) w_i=\frac{1}{\prod_{0\le j\le k,j\not=i}{(x_i-x_j)}} wi=0jk,j=i(xixj)1,则 L ( x ) = l ( x ) ∑ i = 0 k y i w i x − x i L(x)=l(x)\sum_{i=0}^k\frac{y_iw_i}{x-x_i} L(x)=l(x)i=0kxxiyiwi
这样 l i ( x ) l_i(x) li(x) 就能以 O ( 1 ) O(1) O(1) 的时间计算出来,总共耗时 O ( n ) O(n) O(n)

x x x连续取值的情况

当所有 x x x 的取值为连续整数时,基本多项式的分母被 x i x_i xi 分成两段。因此我们与处理好阶乘的乘法逆元,就可以以 O ( n ) O(n) O(n) 的复杂度计算分母了。
上述两种方法结合起来,就能把算法的时间复杂度从 O ( n 2 ) O(n^2) O(n2) 变为 O ( n ) O(n) O(n)

STEP4 经典例题

The Sum of the k-th Powers(CF622F)

There are well-known formulas: ∑ i = 1 n i = 1 + 2 + ⋯ + n = n ( n + 1 ) 2 \sum_{i=1}^ni=1+2+\dots+n=\frac{n(n+1)}{2} i=1ni=1+2++n=2n(n+1), ∑ i = 1 n i 2 = 1 2 + 2 2 + ⋯ + n 2 = n ( 2 n + 1 ) ( n + 1 ) 6 \sum_{i=1}^ni^2=1^2+2^2+\dots+n^2=\frac{n(2n+1)(n+1)}{6} i=1ni2=12+22++n2=6n(2n+1)(n+1) , ∑ i = 1 n i 3 = 1 3 + 2 3 + … n 3 = ( n ( n + 1 ) 2 ) 2 \sum_{i=1}^ni^3=1^3+2^3+\dots_n^3=(\frac{n(n+1)}{2})^2 i=1ni3=13+23+n3=(2n(n+1))2 . Also mathematicians found similar formulas for higher degrees.
Find the value of the sum ∑ i = 1 n i k = 1 k + 2 k + ⋯ + n k \sum_{i=1}^ni^k=1^k+2^k+\dots+n^k i=1nik=1k+2k++nk modulo 1 0 9   +   7 10^9 + 7 109+7 (so you should find the remainder after dividing the answer by the value 1 0 9   +   7 10^9 + 7 109+7).
题目大意: ∑ i = 1 n i k ( m o d 1 0 9 + 7 ) \sum_{i=1}^ni^k(mod 10^9+7) i=1nik(mod109+7),其中 1 ≤ n ≤ 1 0 9 1\le n\le 10^9 1n109 0 ≤ k ≤ 1 0 6 0\le k\le 10^6 0k106

F ( x ) = ∑ i = 1 x i k F(x)=\sum_{i=1}^xi^k F(x)=i=1xik 跟据题目的提示,这是一个 k + 1 k+1 k+1 次多项式函数。因此只要计算 F ( 1 ) , F ( 2 ) , … , F ( k + 2 ) F(1),F(2),\dots,F(k+2) F(1),F(2),,F(k+2),再用拉格朗日插值法就能求得 F ( n ) F(n) F(n)。注意到选取的 x x x 为连续整数,可以用到上面提到的两个优化,总时间复杂度为 O ( k ) O(k) O(k)

#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
int n,k;
ll y[1000005],jc[1000005],l=1ll,w,ans,mod=1e9+7;
ll qp(ll a,int b){
	if(b==0) return 1ll;
	if(b==1) return a;
	ll res=qp(a,b>>1);
	(res*=res)%=mod;
	if(b&1) (res*=a)%=mod;
	return res;
}
ll sub(ll a,ll b){
	return (a>b)?(a-b):(a-b+mod);
}
ll inv(ll a){
	return qp(a,mod-2);
}
void lagrange(){
	for(int i=1;i<=k+2;i+=1){
		(w*=sub(n,i))%=mod;
	}
	jc[0]=1ll;
	for(int i=1;i<=k+2;i+=1){
		jc[i]=(jc[i-1]*i*1ll)%mod;  //阶乘预处理
	}
	for(int i=1;i<=k+2;i+=1){  //拉格朗日插值
		l=(jc[i-1]*jc[k+2-i])%mod;
		if((k+2-i)&1) l=sub(0,w);
		l=inv(l);
		(l*=inv(n-i))%=mod;
		(l*=w)%=mod;
		(l*=y[i])%=mod;
		(ans+=l)%=mod;
	}
	return;
}
int main(){
	scanf("%d%d",&n,&k);
	for(int i=1;i<=k+2;i+=1){
		y[i]=qp(i*1ll,k);
		(y[i]+=y[i-1])%=mod;  //计算y
	}
	if(n<=k+2){
		ans=y[n];  //这种情况直接可以获得答案
	}
	else lagrange();
	printf("%lld\n",ans);
	return 0;
}

谢谢观看,记得点赞

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值