拉格朗日插值
例题 Luogu P4781【模板】拉格朗日插值
给出
个点
,将过这
个点的最多
次的多项式记为
,求
的值。
方法 1:差分法
差分法适用于
的情况。
如,用差分法求某三次多项式
的多项式形式,已知
至
的值分别为
。

第一行为
的连续的前
项;之后的每一行为之前一行中对应的相邻两项之差。观察到,如果这样操作的次数足够多(前提是
为多项式),最终总会返回一个定值,可以利用这个定值求出
的每一项的系数,然后即可将
代入多项式中求解。上例中可求出
。时间复杂度为
。这种方法对给出的点的限制性较强。
方法 2:待定系数法
设
将每个
代入
,有
,这样就可以得到一个由
条
元一次方程所组成的方程组,然后使用 高斯消元 解该方程组求出每一项
,即确定了
的表达式。
如果您不知道什么是高斯消元,请看 高斯消元。
时间复杂度
,对给出点的坐标无要求。
方法 3:拉格朗日插值法
在 多项式部分简介 里我们已经定义了多项式除法。
那么我们会有:

这是显然的,因为
,这个式子显然有
这个因式,所以得证。
这样我们就可以列一个关于
的多项式线性同余方程组:

我们根据中国剩余定理,有:

则
模
意义下的逆元就是:

所以就有:

所以在模意义下
就是唯一的,即:

这就是拉格朗日插值的表达式。
如果要将每一项的系数都算出来,时间复杂度仍为
,但是本题中只用求出
的值,所以在计算上式的过程中直接将
代入即可。

本题中,还需要求解逆元。如果先分别计算出分子和分母,再将分子乘进分母的逆元,累加进最后的答案,时间复杂度的瓶颈就不会在求逆元上,时间复杂度为
。
通常意义下拉格朗日插值的一种推导
由于要求构造一个函数
过点
。首先设第
个点在
轴上的投影为
。
考虑构造
个函数
,使得对于第
个函数
,其图像过
,则可知题目所求的函数
。
那么可以设
,将点
代入可以知道
,所以
。
那么我们就可以从另一个角度推导出通常意义下(而非模意义下)拉格朗日插值的式子为:
。
代码实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31 | #include <cstdio>
const int maxn = 2010;
long long mod = 998244353;
long long n, k, x[maxn], y[maxn], ans, s1, s2;
long long powmod(long long x, long long n) { // 快速幂
long long ret = (long long)1;
while (n) {
if (n % 2 == 1) ret = ret * x % mod;
x = x * x % mod;
n /= 2;
}
return ret;
}
long long inv(long long x) { return powmod(x, mod - 2); } // 求逆元
int main() {
scanf("%lld%lld", &n, &k);
for (int i = 1; i <= n; i++) scanf("%lld%lld", x + i, y + i);
for (int i = 1; i <= n; i++) {
s1 = y[i] % mod;
s2 = (long long)1;
for (int j = 1; j <= n; j++)
if (i != j) s1 = s1 * (k - x[j]) % mod, s2 = s2 * (x[i] - x[j]) % mod;
ans += s1 * inv(s2) % mod;
}
printf("%lld\n", (ans % mod + mod) % mod);
return 0;
}
|
横坐标是连续整数的拉格朗日插值
如果已知点的横坐标是连续整数,我们可以做到
插值。
设要求
次多项式为
,我们已知
(
),考虑代入上面的插值公式:

后面的累乘可以分子分母分别考虑,不难得到分子为:

分母的
累乘可以拆成两段阶乘来算:

于是横坐标为
的插值公式:

预处理
前后缀积、阶乘阶乘逆,然后代入这个式子,复杂度为
。
例题 CF622F The Sum of the k-th Powers
给出
,求
对
取模的值。
本题中,答案是一个
次多项式,因此我们可以线性筛出
的值然后进行
插值。
代码实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52 | // By: [email protected]_er(122461)
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 5, mod = 1e9 + 7;
int n, k, tab[N], p[N], pcnt, f[N], pre[N], suf[N], fac[N], inv[N], ans;
int qpow(int x, int y) {
int ans = 1;
for (; y; y >>= 1, x = 1LL * x * x % mod)
if (y & 1) ans = 1LL * ans * x % mod;
return ans;
}
void sieve(int lim) {
f[1] = 1;
for (int i = 2; i <= lim; i++) {
if (!tab[i]) {
p[++pcnt] = i;
f[i] = qpow(i, k);
}
for (int j = 1; j <= pcnt && 1LL * i * p[j] <= lim; j++) {
tab[i * p[j]] = 1;
f[i * p[j]] = 1LL * f[i] * f[p[j]] % mod;
if (!(i % p[j])) break;
}
}
for (int i = 2; i <= lim; i++) f[i] = (f[i - 1] + f[i]) % mod;
}
int main() {
scanf("%d%d", &n, &k);
sieve(k + 2);
if (n <= k + 2) return printf("%d\n", f[n]) & 0;
pre[0] = suf[k + 3] = 1;
for (int i = 1; i <= k + 2; i++) pre[i] = 1LL * pre[i - 1] * (n - i) % mod;
for (int i = k + 2; i >= 1; i--) suf[i] = 1LL * suf[i + 1] * (n - i) % mod;
fac[0] = inv[0] = fac[1] = inv[1] = 1;
for (int i = 2; i <= k + 2; i++) {
fac[i] = 1LL * fac[i - 1] * i % mod;
inv[i] = 1LL * (mod - mod / i) * inv[mod % i] % mod;
}
for (int i = 2; i <= k + 2; i++) inv[i] = 1LL * inv[i - 1] * inv[i] % mod;
for (int i = 1; i <= k + 2; i++) {
int P = 1LL * pre[i - 1] * suf[i + 1] % mod;
int Q = 1LL * inv[i - 1] * inv[k + 2 - i] % mod;
int mul = ((k + 2 - i) & 1) ? -1 : 1;
ans = (ans + 1LL * (Q * mul + mod) % mod * P % mod * f[i] % mod) % mod;
}
printf("%d\n", ans);
return 0;
}
|
本页面最近更新:2022/9/22 16:20:15,更新历史
发现错误?想一起完善? 在 GitHub 上编辑此页!
本页面贡献者:Ir1d, Marcythm, x4Cx58x54, YanWQ-monad, AtomAlpaca, billchenchina, Chrogeek, Early0v0, EndlessCheng, Enter-tainer, Ghastlcon, Henry-ZHR, hsfzLZH1, kenlig, Peanut-Tang, qwqAutomaton, qz-cqy, rui_er, StudyingFather, swift-zym, Tiphereth-A, TrisolarisHD, Xeonacid
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供,附加条款亦可能应用