拉格朗日插值

例题 Luogu P4781【模板】拉格朗日插值

给出 n 个点 P_i(x_i,y_i) ,将过这 n 个点的最多 n-1 次的多项式记为 f(x) ,求 f(k) 的值。

方法 1:差分法

差分法适用于 x_i=i 的情况。

如,用差分法求某三次多项式 f(x)=\sum_{i=0}^{3} a_ix^i 的多项式形式,已知 f(1) f(6) 的值分别为 1, 5, 14, 30, 55, 91

\begin{array}{cccccccccccc} 1 & & 5 & & 14 & & 30 & & 55 & & 91 & \\ & 4 & & 9 & & 16 & & 25 & & 36 & \\ & & 5 & & 7 & & 9 & & 11 & \\ & & & 2 & & 2 & & 2 & \\ \end{array}

第一行为 f(x) 的连续的前 n 项;之后的每一行为之前一行中对应的相邻两项之差。观察到,如果这样操作的次数足够多(前提是 f(x) 为多项式),最终总会返回一个定值,可以利用这个定值求出 f(x) 的每一项的系数,然后即可将 k 代入多项式中求解。上例中可求出 f(x)=\frac 1 3 x^3+\frac 1 2 x^2+\frac 1 6 x 。时间复杂度为 O(n^2) 。这种方法对给出的点的限制性较强。

方法 2:待定系数法

f(x)=\sum_{i=0}^{n-1} a_ix^i 将每个 x_i 代入 f(x) ,有 f(x_i)=y_i ,这样就可以得到一个由 n n 元一次方程所组成的方程组,然后使用 高斯消元 解该方程组求出每一项 a_i ,即确定了 f(x) 的表达式。

如果您不知道什么是高斯消元,请看 高斯消元

时间复杂度 O(n^3) ,对给出点的坐标无要求。

方法 3:拉格朗日插值法

多项式部分简介 里我们已经定义了多项式除法。

那么我们会有:

f(x)\equiv f(a)\pmod{(x-a)}

这是显然的,因为 f(x)-f(a)=(a_0-a_0)+a_1(x^1-a^1)+a_1(x^2-a^2)+\cdots +a_n(x^n-a^n) ,这个式子显然有 (x-a) 这个因式,所以得证。

这样我们就可以列一个关于 f(x) 的多项式线性同余方程组:

\begin{cases} f(x)\equiv y_1\pmod{(x-x_1)}\\ f(x)\equiv y_2\pmod{(x-x_2)}\\ \cdots\\ f(x)\equiv y_n\pmod{(x-x_n)} \end{cases}

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

M=\prod_{i=1}^n{(x-x_i)},m_i=\dfrac M{x-x_i}=\prod_{j\ne i}{(x-x_j)}

m_i (x-x_i) 意义下的逆元就是:

m_i^{-1}=\prod_{j\ne i}{\dfrac 1{x_i-x_j}}

所以就有:

f(x)\equiv\sum_{i=1}^n{y_im_im_i^{-1}}\equiv\sum_{i=1}^n{y_i\prod_{j\ne i}{\dfrac {x-x_j}{x_i-x_j}}}\pmod M

所以在模意义下 f(x) 就是唯一的,即:

f(x)=\sum_{i=1}^n{y_i\prod_{j\ne i}{\dfrac {x-x_j}{x_i-x_j}}}

这就是拉格朗日插值的表达式。

如果要将每一项的系数都算出来,时间复杂度仍为 O(n^2) ,但是本题中只用求出 f(k) 的值,所以在计算上式的过程中直接将 k 代入即可。

f(k)=\sum_{i=1}^{n} y_i\prod_{j\neq i }\frac{k-x_j}{x_i-x_j}

本题中,还需要求解逆元。如果先分别计算出分子和分母,再将分子乘进分母的逆元,累加进最后的答案,时间复杂度的瓶颈就不会在求逆元上,时间复杂度为 O(n^2)

通常意义下拉格朗日插值的一种推导

由于要求构造一个函数 f(x) 过点 P_1(x_1, y_1), P_2(x_2,y_2),\cdots,P_n(x_n,y_n) 。首先设第 i 个点在 x 轴上的投影为 P_i^{\prime}(x_i,0)

考虑构造 n 个函数 f_1(x), f_2(x), \cdots, f_n(x) ,使得对于第 i 个函数 f_i(x) ,其图像过 \begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases} ,则可知题目所求的函数 f(x)=\sum\limits_{i=1}^nf_i(x)

那么可以设 f_i(x)=a\cdot\prod_{j\neq i}(x-x_j) ,将点 P_i(x_i,y_i) 代入可以知道 a=\dfrac{y_i}{\prod_{j\neq i} (x_i-x_j)} ,所以

f_i(x)=y_i\cdot\dfrac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}=y_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}

那么我们就可以从另一个角度推导出通常意义下(而非模意义下)拉格朗日插值的式子为:

f(x)=\sum_{i=1}^ny_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}

代码实现

 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;
}

横坐标是连续整数的拉格朗日插值

如果已知点的横坐标是连续整数,我们可以做到 O(n) 插值。

设要求 n 次多项式为 f(x) ,我们已知 f(1),\cdots,f(n+1) 1\le i\le n+1 ),考虑代入上面的插值公式:

\begin{aligned} f(x)&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-x_j}{x_i-x_j}\\ &=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-j}{i-j} \end{aligned}

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

\dfrac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i}

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

(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!

于是横坐标为 1,\cdots,n+1 的插值公式:

f(x)=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!}

预处理 (x-i) 前后缀积、阶乘阶乘逆,然后代入这个式子,复杂度为 O(n)

例题 CF622F The Sum of the k-th Powers

给出 n,k ,求 \sum\limits_{i=1}^ni^k 10^9+7 取模的值。

本题中,答案是一个 k+1 次多项式,因此我们可以线性筛出 1^i,\cdots,(k+2)^i 的值然后进行 O(n) 插值。

代码实现
 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: Luogu@rui_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;
}

评论