函数求值 / 【模板】拉格朗日插值
题目链接:ybt金牌导航8-3-1 / luogu P4781
题目大意
给定一个 n 项式上的 n 个点的坐标,然后问你这个多项式 y=f(k) 的 f(k) 值。
思路1
首先,我们考虑拉格朗日插值的朴素法。
那它是怎样的呢?
我们考虑构造
n
n
n 个函数,分别对应
n
n
n 个点,这些函数只在
x
x
x 值为对应点的
x
x
x 值时才为
1
1
1,如果在不是对应点的点,值就为
0
0
0。
那怎么构造出这个函数呢?
可以这样,
k
k
k 是这个函数的
x
x
x 值,
x
,
y
x,y
x,y 是对应点的坐标。
w
(
i
)
=
∏
i
≠
j
k
−
x
j
x
i
−
x
j
w(i)=\prod\limits_{i\neq j}^{}\dfrac{k-x_j}{x_i-x_j}
w(i)=i=j∏xi−xjk−xj
因为如果
k
k
k 是
x
i
x_i
xi,那你会发现,分数上下都是一样,那就都是
1
1
1,乘起来也是
1
1
1。
那如果不是
x
i
x_i
xi,而且是在另一个点上,那它必然会存在一个
j
j
j,使得它的
x
x
x 坐标跟
j
j
j 的一样,那上面就变成了
0
0
0,那这个分数就变成了
0
0
0。有一个乘起来的变成了
0
0
0,那全部也会变成
0
0
0。
那你看如何还原出那个多项式。
这里给出还原方法,其实很好理解:
f
(
k
)
=
y
1
×
w
(
1
)
+
y
2
×
w
(
2
)
+
.
.
.
+
y
n
×
w
(
n
)
f(k)=y_1\times w(1)+y_2\times w(2)+...+y_n\times w(n)
f(k)=y1×w(1)+y2×w(2)+...+yn×w(n)
没错,就是每一个 w w w 函数都是觉得一个点的值为 1 1 1,然后乘上 y y y 值就是那个点的大小。这个时候其它点的 w w w 函数都是 0 0 0,那整个函数的值就是 y y y 值了。
然后就是两个枚举就好了。
f
(
k
)
=
∑
i
=
1
n
(
y
i
∏
i
≠
j
k
−
x
j
x
i
−
x
j
)
f(k)=\sum\limits_{i=1}^{n}(y_i\prod\limits_{i\neq j}\dfrac{k-x_j}{x_i-x_j})
f(k)=i=1∑n(yii=j∏xi−xjk−xj)
如果要还原出多项式,就要 n 3 n^3 n3,不是特别优秀。
代码1
#include<cstdio>
#define ll long long
#define mo 998244353
using namespace std;
int n, k;
ll x[2001], y[2001], ans, sum;
ll ksm(ll x, ll y) {//求逆元的快速幂
ll re = 1;
while (y) {
if (y & 1ll) re = (re * x) % mo;
x = (x * x) % mo;
y >>= 1;
}
return re;
}
int main() {
scanf("%d %d", &n, &k);
for (int i = 1; i <= n; i++) {
scanf("%lld %lld", &x[i], &y[i]);
}
for (int i = 1; i <= n; i++) {
sum = y[i];//乘上的 y[i]
for (int j = 1; j <= n; j++)
if (i != j) {
ll now = (k - x[j]) * ksm(x[i] - x[j], mo - 2) % mo;//最里层的
if (now < 0) now = ((now % mo) + mo) % mo;//可能值会小于0,要把它转回成正的
sum = (sum * now) % mo;//累乘
}
ans = (ans + sum) % mo;//累加
}
printf("%lld", ans);
return 0;
}
思路2
我们考虑弄一下上面的公式。
考虑到那个累乘的
i
≠
j
i\neq j
i=j 很烦,那我们这样,把累乘那里乘上
k
−
x
i
k-x_i
k−xi,然后你的分子就可以整个脱离出来。
然后在前面乘的
y
i
y_i
yi 再除回去。
那就会变成这个:
f
(
k
)
=
∏
i
=
1
n
(
k
−
x
i
)
∑
i
=
1
n
(
y
i
k
−
x
i
∏
i
≠
j
1
x
i
−
x
j
)
f(k)=\prod\limits_{i=1}^{n}(k-x_i)\sum\limits_{i=1}^{n}(\dfrac{y_i}{k-x_i}\prod\limits_{i\neq j}\dfrac{1}{x_i-x_j})
f(k)=i=1∏n(k−xi)i=1∑n(k−xiyii=j∏xi−xj1)
那你如果要求多项式,你可以枚举
i
i
i,然后把
1
k
−
x
i
\dfrac{1}{k-x_i}
k−xi1 拿到外面,然后每次就可以用总共的减去不需要的,就可以
O
(
n
)
O(n)
O(n) 求前面的累乘和你提出来的。
然后其他的也很好求,就可以
O
(
n
2
)
O(n^2)
O(n2) 求。
然后如果要多一个插入的点,你就把最里面的累乘多除以 x i − x n + 1 x_i-x_{n+1} xi−xn+1,然后就可以 O ( n ) O(n) O(n) 求累加的。然后左边的累乘简单一搞就好,那就可以 O ( n ) O(n) O(n) 插入。
而这个东西,就叫做重心拉格朗日插值法。
代码2
#include<cstdio>
#define ll long long
#define mo 998244353
using namespace std;
int n, k;
ll x[2001], y[2001], ans1, sum;
ll ans2;
ll ksm(ll x, ll y) {//快速幂求逆元
if (x < 0) x = (x % mo + mo) % mo;
ll re = 1;
while (y) {
if (y & 1ll) re = (re * x) % mo;
x = (x * x) % mo;
y >>= 1;
}
return re;
}
int main() {
scanf("%d %d", &n, &k);
for (int i = 1; i <= n; i++) {
scanf("%lld %lld", &x[i], &y[i]);
}
ans1 = 1;
for (int i = 1; i <= n; i++) {//前面的累乘
ll now = k - x[i];
if (now < 0) now = (now % mo + mo) % mo;
ans1 = (ans1 * now) % mo;
}
for (int i = 1; i <= n; i++) {//按着公式算
ll sum = (y[i] * ksm(k - x[i], mo - 2)) % mo;
for (int j = 1; j <= n; j++)
if (i != j)
sum = (sum * ksm(x[i] - x[j], mo - 2)) % mo;
ans2 = (ans2 + sum) % mo;
}
printf("%lld", (ans1 * ans2) % mo);
return 0;
}
一些小小的补充
这里是求一个点的,如果你要求多个点的时候,你可以提前预处理:
t
i
=
∏
i
≠
j
1
x
i
−
x
j
t_i=\prod\limits_{i\neq j}\dfrac{1}{x_i-x_j}
ti=i=j∏xi−xj1(
O
(
n
2
)
O(n^2)
O(n2) 预处理)
然后式子就是:
f
(
k
)
=
∏
i
=
1
n
(
k
−
x
i
)
∑
i
=
1
n
(
y
i
t
i
k
−
x
i
)
f(k)=\prod\limits_{i=1}^{n}(k-x_i)\sum\limits_{i=1}^{n}(\dfrac{y_it_i}{k-x_i})
f(k)=i=1∏n(k−xi)i=1∑n(k−xiyiti)
然后就是
O
(
n
)
O(n)
O(n) 算每个点了。