拉格朗日插值法

拉格朗日插值法

简介

定义：给定 \(n\) 个点 \((x_i,y_i)\)，求一个次数不超过 \(n-1\) 的多项式函数 \(f(x)\) 使得 \(f(x_i)=y_i\)。要求 \(x_i\) 两两互不相同。

下文注意区分函数变量 \(x\) 与给定点坐标 \(x_i\)。

拉格朗日插值基函数：定义 \(f_i(x)=\prod_{i\not=j}\frac{x-x_j}{x_i-x_j}\)。这样定义的好处是，当代入 \(x_j\not=x_i\) 时，\(f_i(x_j)=0\)；当代入 \(x_i\) 时，\(f_i(x_i)=y_i\)。

那么要求即为函数 \(f(x)=\sum_{i=1}^{n}f_i(x)\)。

有定理：\(n\) 个横坐标互不相同的点能确定一个次数小于等于 \(n-1\) 的函数。存在性由拉格朗日插值法给出。唯一性可以考虑反证，如果存在两个满足条件的函数，那么将两函数相减，得到的函数次数不超过 \(n-1\)，但至少有 \(n\) 个零点，与代数基本定理矛盾。

插值求值

给定 \(n\) 个点 \((x_i,y_i)\)，求在 \(x=k\) 处的取值。只需要将上面的式子中的变量 \(x\) 直接代入即可。时间复杂度为 \(O(n^2)\)。

Code

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
using ldb = long double;

#define file(a) freopen(a ".in", "r", stdin), freopen(a ".out", "w", stdout)
#define unsyncio() cin.tie(nullptr)->sync_with_stdio(false)
#define eb emplace_back
#define mpi make_pair

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
template <typename T = int>
T rad(T l, T r) {
  return uniform_int_distribution<T>(l, r)(rng);
}

const int N = 2e3 + 5, mod = 998244353;

inline i64 qpow(i64 a, i64 b) {
  i64 res = 1;
  for (; b; b >>= 1, a = a * a % mod)
    if (b & 1)
      res = res * a % mod;
  return res;
}
inline i64 inv(i64 a) {
  return qpow(a, mod - 2);
}
int x[N], y[N];

int main() {
  int n, xx;
  cin >> n >> xx;
  for (int i = 1; i <= n; ++i) {
    cin >> x[i] >> y[i];
  }
  i64 ans = 0;
  for (int i = 1; i <= n; ++i) {
    i64 v = 1;
    for (int j = 1; j <= n; ++j) {
      if (i != j) {
        v = v * (xx + mod - x[j]) % mod;
        v = v * inv(x[i] + mod - x[j]) % mod;
      }
    }
    ans = (ans + v * y[i] % mod) % mod;
  }
  cout << ans << '\n';
  return 0;
}

求系数

回忆一下刚刚的式子：

\[f(x)=\sum_{i=1}^{n}y_i\prod_{j\not=i}\frac{x-x_j}{x_i-x_j} \]

将分子分母同时乘上 \(x-x_i\)，分子全部变为 \(g(x)=\prod_{i=1}^{n}(x-x_i)\)，即

\[f(x)=\sum_{i=1}^{n}y_i\cdot\frac{g(x)}{x-x_i}\prod_{j\not=i}\frac{1}{x_i-x_j} \]

这启示我们先用 \(O(n^2)\) 暴力乘出 \(g(x)\) 的系数表达——这不困难。然后再对每一个 \(i\) 求出其除以 \(x-x_i\) 的系数表达——这也不困难。最后乘上常数，再加起来就可以了。

Code

// 拉格朗日插值法（求系数）

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
using ldb = long double;

#define file(a) freopen(a ".in", "r", stdin), freopen(a ".out", "w", stdout)
#define unsyncio() cin.tie(nullptr)->sync_with_stdio(false)
#define eb emplace_back
#define mpi make_pair

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
template <typename T = int>
T rad(T l, T r) {
  return uniform_int_distribution<T>(l, r)(rng);
}

const int N = 2e3 + 5, mod = 998244353;
int x[N], y[N];

i64 qpow(i64 a, i64 b) {
  i64 res = 1;
  for (; b; b >>= 1, a = a * a % mod)
    if (b & 1)
      res = res * a % mod;
  return res;
}

int main() {
  unsyncio();
  int n, k;
  cin >> n >> k;
  for (int i = 1; i <= n; ++i) {
    cin >> x[i] >> y[i];
  }
  static int bas[N];
  bas[0] = 1;
  for (int i = 1; i <= n; ++i) {
    for (int j = i; j >= 1; --j) {
      bas[j] = (bas[j - 1] + mod - 1ull * bas[j] * x[i] % mod) % mod;
    }
    bas[0] = (mod - 1ull * bas[0] * x[i] % mod) % mod;
  }
  static int ans[N];
  for (int i = 1; i <= n; ++i) {
    static int tmp[N];
    for (int j = n, res = 0; j >= 0; --j) {
      tmp[j] = res;
      if (j)
        res = (bas[j] + 1ull * res * x[i] % mod) % mod;
    }
    int val = 1;
    for (int j = 1; j <= n; ++j)
      if (i != j)
        val = 1ull * val * (x[i] + mod - x[j]) % mod;
    val = 1ull * qpow(val, mod - 2) * y[i] % mod;
    for (int j = 0; j < n; ++j) {
      (ans[j] += 1ull * tmp[j] * val % mod) %= mod;
    }
  }
  int anss = 0;
  for (int i = 0, p = 1; i < n; ++i) {
    (anss += 1ull * ans[i] * p % mod) %= mod;
    p = 1ull * p * k % mod;
  }
  cout << anss << '\n';
  return 0;
}

重心拉格朗日插值法

有点类似于上面提到的内容。

\[f(x)=\sum_{i=1}^{n}y_i\cdot\frac{g(x)}{x-x_i}\prod_{j\not=i}\frac{1}{x_i-x_j} \]

令这个式子右边 \(w_i=\prod_{j\not=i}\frac{1}{x_i-x_j}\)，这个 \(w_i\) 被称为重心权。

有人发现这个 \(w_i\) 在插入一个新的点的时候可以 \(O(n\log n)\) 修改，这就实现了 \(O(n\log n)\) 插入。

单次求值应该可以做到 \(O(n\log n)\)：先求 \(g(k)\)，再每一个除以 \(k-x_i\)。

但是我在实现的时候发现，求值得到 \(g(k)\) 时可能出现 \(k=x_i\) 的情况，我特判掉直接输出 \(y_i\) 了，不知道有没有更好的解决方案。

Code

// 重心拉格朗日插值法

#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
using u32 = unsigned int;
using u64 = unsigned long long;
using ldb = long double;

#define file(a) freopen(a ".in", "r", stdin), freopen(a ".out", "w", stdout)
#define unsyncio() cin.tie(nullptr)->sync_with_stdio(false)
#define eb emplace_back
#define mpi make_pair

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
template <typename T = int>
T rad(T l, T r) {
  return uniform_int_distribution<T>(l, r)(rng);
}

const int N = 2e3 + 5, mod = 998244353;
int w[N], x[N], y[N];

i64 qpow(i64 a, i64 b) {
  i64 res = 1;
  for (; b; b >>= 1, a = a * a % mod)
    if (b & 1)
      res = res * a % mod;
  return res;
}

int main() {
  unsyncio();
  int n, k;
  cin >> n >> k;
  for (int i = 1; i <= n; ++i) {
    cin >> x[i] >> y[i];
    w[i] = 1;
    for (int j = 1; j < i; ++j) {
      i64 tmp = qpow(x[j] + mod - x[i], mod - 2);
      w[j] = 1ull * w[j] * tmp % mod;
      w[i] = 1ull * w[i] * (mod - tmp) % mod;
    }
  }
  i64 g = 1;
  for (int i = 1; i <= n; ++i) {
    if (k == x[i]) {
      cout << y[i] << '\n';
      return 0;
    }
    g = 1ull * g * (k + mod - x[i]) % mod;
  }
  i64 ans = 0;
  for (int i = 1; i <= n; ++i) {
    (ans += 1ull * g * qpow(k + mod - x[i], mod - 2) % mod * w[i] % mod * y[i] % mod) %= mod;
  }
  cout << ans << '\n';
  return 0;
}