拉格朗日插值法
无关内容
- 无关前置:今天学习
L
a
g
r
a
n
g
e
\rm Lagrange
Lagrange,然后我一直疑惑为什么音译过来是拉格朗日(不应该是拉格阮籍吗?QAQ),最后在同学的帮助下,在知乎上找到了原因:原来真相居然是,法语,emmmm
进入正题
现在给你一个 n n n次多项式的点值表示法(说简单点就是给你一个多项式函数图像上的 n + 1 n+1 n+1个点对 ( x i , y i ) (x_i,y_i) (xi,yi)),请你求出当 x = k x=k x=k时的函数值(也就是将 k k k带入多项式求值)。
-
暴力想法:我们有了 n + 1 n+1 n+1个点对,那么将其分别带入可以得到 n + 1 n+1 n+1个 n n n元一次(本来是一元 n n n次方程,但是将 x r x^r xr看作不同的未知数就变成了 n + 1 n+1 n+1元一次方程),然后高斯消元即可,复杂度 O ( n 3 ) O(n^3) O(n3)。
-
进一步思考:但是对于 n n n达到几千的情况, O ( n 3 ) O(n^3) O(n3)显然会超时,那么有没有更快的方法呢?
-
我们可以用斜率来拟合原函数,我们可以通过先枚举一个点,然后求给定的 x = k x=k x=k点与其他点的 x i x_i xi组成的两条直线的斜率之比:
假设插入的 k k k求出点对为 ( x , y ) (x,y) (x,y),现在枚举的 i , j ( i ≠ j ) i,j(i\neq j) i,j(i̸=j) , 那么斜率之比为 x − x j y − y j x i − x j y i − y j \frac{\frac{x-x_j}{y-y_j}}{\frac{x_i-x_j}{y_i-y_j}} yi−yjxi−xjy−yjx−xj,整理得 x − x j x i − x j × y i − y j y − y j \frac{x-x_j}{x_i-x_j}\times \frac{y_i-y_j}{y-y_j} xi−xjx−xj×y−yjyi−yj,那么比例就大概为 x − x j x i − x j \frac{x-x_j}{x_i-x_j} xi−xjx−xj(因为y并不知道)
然后用直线的斜率之比的乘积来当做新的斜率之比,然后我们就拟合出了一条 ∏ i ≠ j ( x − x j x i − x j ) \prod\limits_{i\neq j}(\frac{x-x_j}{x_i-x_j}) i̸=j∏(xi−xjx−xj)斜率的直线,那么该点对于给定的 x = k x=k x=k的贡献(也就是对于的 y y y值)可以由比例得知为 y i ∏ i ≠ j ( x − x j x i − x j ) y_i\prod\limits_{i\neq j}(\frac{x-x_j}{x_i-x_j}) yii̸=j∏(xi−xjx−xj),那么将 n + 1 n+1 n+1个点的所有贡献累加起来便得知了当 x = k x=k x=k的时候 y = ∑ i = 0 n y i ∏ i ≠ j ( x − x j x i − x j ) y=\sum\limits_{i=0}^{n}y_i\prod\limits_{i\neq j}(\frac{x-x_j}{x_i-x_j}) y=i=0∑nyii̸=j∏(xi−xjx−xj)( i ≠ j i\neq j i̸=j是因为自己和自己一个点是不能求斜率的)
所以,我们正式推出拉格朗日插值公式:
∑
i
=
0
n
y
i
∏
i
≠
j
(
x
−
x
j
x
i
−
x
j
)
\sum\limits_{i=0}^{n}y_i\prod\limits_{i\neq j}(\frac{x-x_j}{x_i-x_j})
i=0∑nyii̸=j∏(xi−xjx−xj)
拉格朗日插值法可能会有一定的误差,因为是拟合出来的。
(证明与推导可能有些不妥或错误的地方,可以提出或者参考其他的blog或者百科)
此公式还有另外一个证明,就是将给出的
x
i
x_i
xi带入该公式,我们会发现,出除了当你枚举的
i
i
i等于带入的
x
i
x_i
xi的
i
i
i时,其他情况都有一个
x
i
−
x
i
=
0
x_i-x_i=0
xi−xi=0,而当相等时,则有
y
i
∏
i
≠
j
(
x
i
−
x
j
x
i
−
x
j
)
y_i\prod\limits_{i\neq j}(\frac{x_i-x_j}{x_i-x_j})
yii̸=j∏(xi−xjxi−xj),我们会发现左边的值一定为
1
1
1(分子分母相同),那么最后插值的答案确实就等于
y
i
y_i
yi,所以是正确的。
推广来说,我们令
β
i
(
x
)
=
∏
i
≠
j
x
−
x
j
x
i
−
x
j
\beta_i(x)=\prod\limits_{i\neq j}\frac{x-x_j}{x_i-x_j}
βi(x)=i̸=j∏xi−xjx−xj,我们就可以得知函数
β
i
(
x
i
)
=
1
\beta_i(x_i)=1
βi(xi)=1,而
β
i
(
x
j
)
=
0
\beta_i(x_j)=0
βi(xj)=0
那么,我们通过公式可以看出,这个方法的复杂度是 O ( n 2 ) O(n^2) O(n2)的,如果涉及取模和逆元,那么使用快速幂,复杂度则为 O ( n l o g v a l + n 2 ) O(n_{log}val+n^2) O(nlogval+n2)还是 O ( n 2 ) O(n^2) O(n2)。
代码实现就非常简单啦!
丑陋代码如下:
点对为x[i],y[i]。
n-1次多项式,求插入k的值
#define ll long long
#define fpow(a,b) pow(a,b)
ll Lagrange(ll *x,ll *y,int n,ll k){
ll s1,s2,ans=0;
for(RG int i=1;i<=n;++i){
s1=1;s2=1;
for(RG int j=1;j<=n;++j){
if(i==j) continue;
s1=(s1*(k-x[j]+Mod)%Mod)%Mod;
s2=(s2*(x[i]-x[j]+Mod)%Mod)%Mod;
}
ans=(ans+s1*fpow(s2,Mod-2)%Mod*y[i]%Mod)%Mod;
}
return ans;
}
- 特殊情况:
有时我们的
x
i
x_i
xi会是连续的取值,如
x
i
x_i
xi为
1
∼
n
1\sim n
1∼n的连续正整数,这时再观察原式子,则变成了:
∑
i
=
0
n
y
i
∏
i
≠
j
x
−
j
i
−
j
\sum\limits_{i=0}^ny_i \prod \limits_{i\neq j} \frac{x-j}{i-j}
i=0∑nyii̸=j∏i−jx−j
于是我们对于分子预处理前缀积和后缀积,对于分母再预处理阶乘(取模还有阶乘的逆元),那么就 O ( n + l o g M a x v a l + n ) O(n+logMaxval+n) O(n+logMaxval+n)相当于 O ( n ) O(n) O(n)求出来了。
此时插入
x
=
k
x=k
x=k:
前缀积:
p
r
e
[
i
]
=
∏
j
=
1
i
(
k
−
j
)
pre[i]=\prod\limits_{j=1}^i(k-j)
pre[i]=j=1∏i(k−j)
后缀积:
s
u
f
[
i
]
=
∏
j
=
i
n
(
k
−
j
)
suf[i]=\prod\limits_{j=i}^n(k-j)
suf[i]=j=i∏n(k−j)
那么,显然分子就等于
p
r
e
[
i
−
1
]
×
s
u
f
[
i
+
1
]
pre[i-1]\times suf[i+1]
pre[i−1]×suf[i+1](刚好取出第
i
i
i个)
阶乘:
f
a
c
[
i
]
=
∏
j
=
1
i
j
fac[i]=\prod\limits_{j=1}^ij
fac[i]=j=1∏ij
逆元可以这样处理:(这里使用费马小定理,模数为质数,你也可以使用其他定理求))
i
n
v
f
a
c
[
n
]
=
p
o
w
(
f
a
c
[
n
]
,
m
o
d
−
2
)
invfac[n]=pow(fac[n],mod-2)
invfac[n]=pow(fac[n],mod−2)
然后线性递推出所有的逆元,
i
n
v
f
a
c
[
i
]
=
i
n
v
f
a
c
[
i
+
1
]
×
(
i
+
1
)
invfac[i]=invfac[i+1]\times (i+1)
invfac[i]=invfac[i+1]×(i+1)
那么,显然分母就等于
i
n
v
f
a
c
[
i
−
1
]
×
i
n
v
f
a
c
[
n
−
i
]
invfac[i-1]\times invfac[n-i]
invfac[i−1]×invfac[n−i]
但是这里分母还有一个正负的问题,其实不难发现,当 n − i n-i n−i为奇数时就是负数,否则为偶数。
那么代码也就显然啦,下面上丑陋代码QAQ
ll fac_inv[M];
ll pre[M],suf[M];
ll SpLagrange(/*ll *x,(x为1~n连续整数)*/ll *y,int n,ll k){
pre[0]=1;for(RG int i=1;i<=n;++i)pre[i]=pre[i-1]*((k-i+Mod)%Mod)%Mod;
suf[n+1]=1;for(RG int i=n;i>=1;--i)suf[i]=suf[i+1]*((k-i+Mod)%Mod)%Mod;
ll fac=1;for(RG int i=2;i<=n;++i)fac=fac*i%Mod;
fac=fpow(fac,Mod-2);fac_inv[n]=fac;
for(RG int i=n-1;i>=1;--i)fac_inv[i]=fac_inv[i+1]*(i+1)%Mod;
for(RG int i=1;i<=n;++i){
ll s1=pre[i-1]*suf[i+1]%Mod;
ll s2=inv_fac[i-1]*inv_fac[n-i]%Mod;
ll flag=((n-i)&1)?-1:1;
ans=(ans+flag*s1*s2%Mod*y[i]%Mod)%Mod;
}
return (ans+Mod)%Mod;
}
- 既然已经有了一些基础知识,那么我们来做几道例题吧:
自然数的幂的前缀和
求下列式子的值:
∑
i
=
0
n
i
k
\sum\limits_{i=0}^ni^k
i=0∑nik
数据范围: 1 ≤ n ≤ 1 0 15 , 0 ≤ k ≤ 1 0 6 1\leq n\leq10^{15},0\leq k\leq 10^6 1≤n≤1015,0≤k≤106, M o d = 998244353 Mod=998244353 Mod=998244353
这是拉格朗日插值法的经典题目 C F 622 F T h e S u m o f t h e k t h P o w e r s \rm CF622F\ The\ Sum\ of\ the\ k^{th}\ Powers CF622F The Sum of the kth Powers,根据大佬说,也可以使用伯努利数 + 任意模数NTT(大佬文章【真正讲解伯努利数的文章】),但是蒟蒻并不会QAQ(要用到生成函数和NTT的知识),所以我们还是老老实实用拉格朗日插值法吧(而且复杂度还要小一些)。
首先,显然复杂度不能有 n n n,(最多也只能 l o g n logn logn),所以我们要从较小的 k k k上分析。
通过经验,我们知道:
当
k
=
0
k=0
k=0时,求值公式为
n
n
n,一个一次多项式。
当
k
=
1
k=1
k=1时,求值公式为
n
×
(
n
+
1
)
2
\frac{n\times (n+1)}{2}
2n×(n+1),一个二次多项式。
当
k
=
2
k=2
k=2时,求值公式为
n
×
(
n
+
1
)
×
(
2
n
+
1
)
6
\frac{n\times (n+1)\times (2n+1)}{6}
6n×(n+1)×(2n+1),一个三次多项式。
当
k
=
3
k=3
k=3时,求值公式为
(
n
×
(
n
+
1
)
2
)
2
\left(\frac{n\times (n+1)}{2}\right)^2
(2n×(n+1))2,一个四次多项式。
当
k
=
4
k=4
k=4时,求值公式为
n
(
n
+
1
)
(
2
n
+
1
)
(
3
n
2
+
3
n
−
1
)
30
\frac{n(n+1)(2n+1)(3n^2+3n-1)}{30}
30n(n+1)(2n+1)(3n2+3n−1),一个五次多项式。
⋯
\cdots
⋯
由此,我们大胆猜测,
∑
i
=
0
n
i
k
\sum_{i=0}^ni^k
∑i=0nik为一个
k
+
1
k+1
k+1次多项式。
证明(方法来自 I m a g i n e \rm Imagine Imagine大佬的blog):
-
我们将两个 k + 1 k+1 k+1次多项式 ( n + 1 ) k + 1 (n+1)^{k+1} (n+1)k+1与 n k + 1 n^{k+1} nk+1相减,结合二项式定理展开公式,得到( ( n m ) \tbinom{n}{m} (mn)表示组合数):
( n + 1 ) k + 1 − n k + 1 = ∑ i = 0 k + 1 ( k + 1 i ) n i − n k + 1 = ∑ i = 0 k ( k + 1 i ) n i (n+1)^{k+1}-n^{k+1}=\sum\limits_{i=0}^{k+1} \tbinom{k+1}{i}n^i-n^{k+1}\\ =\sum\limits_{i=0}^k\tbinom{k+1}{i}n^i (n+1)k+1−nk+1=i=0∑k+1(ik+1)ni−nk+1=i=0∑k(ik+1)ni -
再将 n k + 1 n^{k+1} nk+1减去 ( n − 1 ) k + 1 (n-1)^{k+1} (n−1)k+1,同理得到:
n k + 1 − ( n − 1 ) k + 1 = ∑ i = 0 k ( k + 1 i ) ( n − 1 ) i n^{k+1}-(n-1)^{k+1}=\sum\limits_{i=0}^k\tbinom{k+1}{i}(n-1)^i nk+1−(n−1)k+1=i=0∑k(ik+1)(n−1)i -
我们令 ζ ( x ) = ∑ i = 0 k ( k + 1 i ) ( n − x ) i \zeta(x)=\sum\limits_{i=0}^k\tbinom{k+1}{i}(n-x)^i ζ(x)=i=0∑k(ik+1)(n−x)i,然后我们求出下列式子:
∑ x = − 1 n ζ ( x ) \sum_{x=-1}^n\zeta(x) x=−1∑nζ(x)
化简得到 ( n + 1 ) k + 1 = ∑ i = 0 k ( k + 1 i ) i k (n+1)^{k+1}=\sum\limits_{i=0}^k\tbinom{k+1}{i}i^k (n+1)k+1=i=0∑k(ik+1)ik,那么显然 ∑ i = 0 n i k \sum_{i=0}^ni^k ∑i=0nik为一个 k + 1 k+1 k+1次多项式。
那么原式是一个
k
+
1
k+1
k+1次的多项式,所以我们取
k
+
2
k+2
k+2个值
x
i
x_i
xi,并计算出对应取值
y
i
=
∑
i
=
1
x
i
i
k
y_i=\sum_{i=1}^{x_i}i^k
yi=∑i=1xiik,显然,我们发现当你的
x
i
x_i
xi取
[
1
,
k
+
2
]
[1,k+2]
[1,k+2]中的正整数值时最简单,此时这个
x
i
x_i
xi也符合前面讲的特殊情况,所以我们可以先预处理求出
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi),然后带入下列式子:
∑
i
=
1
n
i
k
=
∑
i
=
1
k
+
2
y
i
∏
j
=
1
,
j
≠
i
k
+
2
n
−
j
i
−
j
\sum\limits_{i=1}^ni^k=\sum\limits_{i=1}^{k+2}y_i\prod\limits_{j=1,j\neq i}^{k+2}\frac{n-j}{i-j}
i=1∑nik=i=1∑k+2yij=1,j̸=i∏k+2i−jn−j
最后的复杂度为
O
(
k
l
o
g
k
)
O(klogk)
O(klogk)的预处理和
O
(
k
)
O(k)
O(k)的插值。
下面提供几道例题:
BZOJ4559 别人的讲解
HDU4059 别人的讲解
HDU6239 别人的讲解
其他的较好的文章 :