拉格朗日插值法
1.拉格朗日插值法的由来
拉格朗日插值法究竟解决的是什么问题喃?
问题:在二维图像上给出 n n n个点,如何得到一个通过这 n n n个点的一条线?
例如:已知以下几个点 ( 0 , 1 ) 、 ( 1 , 0 ) 、 ( 2 , 1 ) (0,1)、(1,0)、(2,1) (0,1)、(1,0)、(2,1),问这条曲线是什么?
显然易见对于这三个点,是可以用多项式画出这三个点的,因此我们可以假设这个二次多项式为
y
=
a
0
+
a
1
x
+
a
2
x
2
y=a_0+a_1x+a_2x^2
y=a0+a1x+a2x2
我们带入这三个点
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
(
x
3
,
y
3
)
(x_1,y_1),(x_2,y_2),(x_3,y_3)
(x1,y1),(x2,y2),(x3,y3)
可以得到以下方程组
{
y
1
=
a
0
+
a
1
x
1
+
a
2
x
1
2
y
2
=
a
0
+
a
1
x
2
+
a
2
x
2
2
y
3
=
a
0
+
a
1
x
3
+
a
2
x
3
2
\begin{cases} y_1=a_0+a_1x_1+a_2x_1^2\\y_2=a_0+a_1x_2+a_2x_2^2\\y_3=a_0+a_1x_3+a_2x_3^2\\ \end{cases}
⎩
⎨
⎧y1=a0+a1x1+a2x12y2=a0+a1x2+a2x22y3=a0+a1x3+a2x32
对于这样一个方程组,我们可以用高斯消元解出系数。
但是拉格朗日是如何思考这个问题的喃?
拉格朗日的思考:
拉格朗日从另外的角度思考了这个问题,在他看来可以通过三根二次曲线相加来达到目标
。那这是怎么的三根二次曲线呢?
这三根曲线只需满足一些条件即可:
曲线
1
:
y
1
f
1
(
x
)
可以保证,在
x
1
点处,取值为
y
1
,
其余两点取值为
0
曲线
2
:
y
2
f
1
(
x
)
可以保证,在
x
2
点处,取值为
y
2
,
其余两点取值为
0
曲线
3
:
y
3
f
1
(
x
)
可以保证,在
x
3
点处,取值为
y
3
,
其余两点取值为
0
曲线1:y_1f_1(x)可以保证,在x_1点处,取值为y_1,其余两点取值为0\\ 曲线2:y_2f_1(x)可以保证,在x_2点处,取值为y_2,其余两点取值为0\\ 曲线3:y_3f_1(x)可以保证,在x_3点处,取值为y_3,其余两点取值为0
曲线1:y1f1(x)可以保证,在x1点处,取值为y1,其余两点取值为0曲线2:y2f1(x)可以保证,在x2点处,取值为y2,其余两点取值为0曲线3:y3f1(x)可以保证,在x3点处,取值为y3,其余两点取值为0
则所求曲线为:
f
(
x
)
=
y
1
f
1
(
x
)
+
y
2
f
2
(
x
)
+
y
3
f
3
(
x
)
f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)
f(x)=y1f1(x)+y2f2(x)+y3f3(x)
而这样一条曲线正好可以过给的三个点
2、拉格朗日插值法的推导
对于过每一个点 ( x i , y i ) (x_i,y_i) (xi,yi)的曲线的函数设为 y j f i ( x j ) y_jf_i(x_j) yjfi(xj)
因此需要满足的条件为:
f
i
(
x
j
)
=
{
1
,
i
=
j
0
,
i
≠
j
f_i(x_j)= \begin{cases} 1,i=j \\ 0,i\neq j \end{cases}
fi(xj)={1,i=j0,i=j
因此我们最后构造出
f
i
(
x
j
)
f_i(x_j)
fi(xj)即可
显而易见,以下这个函数满足条件
f
i
(
x
)
=
(
x
−
x
1
)
(
x
−
x
2
)
.
.
.
(
x
−
x
n
)
(
x
i
−
x
1
)
(
x
i
−
x
2
)
.
.
.
(
x
i
−
x
n
)
(
i
≠
j
)
f_i(x)=\frac{(x-x_1)(x-x_2)...(x-x_n)}{(x_i-x_1)(x_i-x_2)...(x_i-x_n)}(i\neq j)
fi(x)=(xi−x1)(xi−x2)...(xi−xn)(x−x1)(x−x2)...(x−xn)(i=j)
因此
f
(
x
)
=
∑
i
=
1
n
y
i
f
i
(
x
)
=
∑
i
=
1
n
y
i
Π
i
≠
j
x
−
x
j
x
i
−
x
j
f(x)=\sum_{i=1}^{n}y_if_i(x)=\sum_{i=1}^{n}y_i\Pi_{i\neq j}\frac{x-x_j}{x_i-x_j}
f(x)=i=1∑nyifi(x)=i=1∑nyiΠi=jxi−xjx−xj
上模板题:https://www.luogu.com.cn/problem/P4781
#include<iostream>
#include<algorithm>
#include<cmath>
#include<string.h>
using namespace std;
#define int long long
const int maxn=1e4+10;
int x[maxn];
int y[maxn];
const int mod =998244353;
int qp(int x,int n){
int ans=1;
while(n){
if(n&1){
ans=ans*x%mod;
}
x=x*x%mod;
n>>=1;
}
return ans;
}
int inv(int x){
return qp(x,mod-2);
}
int lagrange(int num,int n)
{
int ans=0;
for(int k=1;k<=n;k++){
int s1=y[k];
int s2=1;
for(int i=1;i<=n;i++){
if(i==k)continue;
s1=s1*(num-x[i])%mod;
s2=s2*(x[k]-x[i])%mod;
}
ans=(ans+s1*inv(s2))%mod;
}
ans=(ans%mod+mod)%mod;
return ans;
}
signed main()
{
int n,k;
cin>>n>>k;
for(int i=1;i<=n;i++){
cin>>x[i]>>y[i];
}
cout<<lagrange(k,n)<<endl;
return 0;
}
当然这样的时间复杂度 O ( n 2 ) O(n^2) O(n2)
在有些情况下不能做到最优,因此对于一些情况我们可以进行优化
3、横坐标为连续的自然数(1,2,3…)时,可以进行优化
由于横坐标为连续的自然数:
故插值公式可以优化为:
f
(
x
)
=
∑
i
=
1
n
+
1
y
i
Π
j
=
1
n
+
1
(
x
−
j
)
(
x
−
j
)
(
−
1
)
n
+
1
−
i
(
i
−
1
)
!
(
n
+
1
−
i
)
!
f(x)=\sum_{i=1}^{n+1}y_i\frac{\Pi_{j=1}^{n+1}(x-j)}{(x-j)(-1)^{n+1-i}(i-1)!(n+1-i)!}
f(x)=i=1∑n+1yi(x−j)(−1)n+1−i(i−1)!(n+1−i)!Πj=1n+1(x−j)
这样一来,可以通过预处理后缀积、前缀积然后代入这个式子,所以这种情况时间复杂度为
O
(
n
)
O(n)
O(n)
例题:求自然数的k次幂的和
求
∑
i
=
1
n
i
k
思路:首先我们通过这个式子可以判断出这个式子最终一定是一个
k
+
1
次的多项式
因此我们只需要
k
+
2
个点
那么接下来如何设计这个函数喃?我们可以把所求看为
f
(
n
)
=
∑
i
=
1
n
i
k
那么每一个点为
(
j
,
∑
i
=
1
j
i
k
)
因此我们预处理出
k
+
2
个点,然后将这
k
+
2
个点代入插值公式即可
求\sum_{i=1}^{n}i^k\\ 思路:首先我们通过这个式子可以判断出这个式子最终一定是一个k+1次的多项式\\ 因此我们只需要k+2个点\\ 那么接下来如何设计这个函数喃?我们可以把所求看为f(n)=\sum_{i=1}^{n}i^k\\ 那么每一个点为(j,\sum_{i=1}^{j}i^k)\\ 因此我们预处理出k+2个点,然后将这k+2个点代入插值公式即可
求i=1∑nik思路:首先我们通过这个式子可以判断出这个式子最终一定是一个k+1次的多项式因此我们只需要k+2个点那么接下来如何设计这个函数喃?我们可以把所求看为f(n)=i=1∑nik那么每一个点为(j,i=1∑jik)因此我们预处理出k+2个点,然后将这k+2个点代入插值公式即可
相同题意的例题:https://codeforces.com/problemset/problem/622/F
代码:
#include <iostream>
#include <algorithm>
#include <cmath>
#include <string.h>
#define int long long
using namespace std;
const int maxn = 1e6 + 10;
const int mod = 1e9 + 7;
int x[maxn];
int y[maxn];
int qp(int a, int b)
{
int ans = 1;
while (b)
{
if (b & 1)
{
ans = ans * a % mod;
}
a = a * a % mod;
b >>= 1;
}
return ans;
}
int inv(int x)
{
return qp(x, mod - 2);
}
int fac[maxn];
void init_fac()
{
fac[0] = 1;
fac[1] = 1;
for (int i = 2; i < maxn - 5; i++)
{
fac[i] = fac[i - 1] * i % mod;
}
}
int back;
int lagrange(int n, int k)
{
int ans = 0;
for (int i = 1; i <= k + 2; i++)
{
int s1 = back * y[i] % mod;
int s2 = ((k + 2 - i) % 2 == 0 ? 1 : -1) * (n - i);
s2 = s2 * fac[i - 1] % mod;
s2 = s2 * fac[k + 2 - i] % mod;
ans = (ans + s1 * inv(s2) % mod) % mod;
}
ans = (ans % mod + mod) % mod;
return ans;
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
int n, k;
cin >> n >> k;
init_fac();
for (int i = 1; i <= k + 2; i++)
{
x[i] = i;
y[i] = (y[i - 1] + qp(i, k)) % mod;
}
if (n <= k+2)
{
cout << y[n] << endl;
}
else
{
back = 1;
for (int j = 1; j <= k + 2; j++)
{
back = back * (n - j) % mod;
}
cout << lagrange(n, k) << endl;
}
return 0;
}