拉格朗日插值

yabnto發表於2024-09-07

CSR：又拉又插的東西（又垃圾，又傻叉）

JCY：什麼你拉插了一晚上？（我學習拉插學了一晚上）

什麼是拉插

給定一些點值，是否可以求出一個函式，使得函式影像穿過這些點，並求出給定的 \(x\) 所對應的 \(y\)。

初步思路

前置

我們想到兩個點肯定可以確定一條直線，而三個點肯定可以確定一條拋物線，所以我們猜測 \(n + 1\) 個點肯定可以確定一條 \(n\) 次的函式

正文：構造

考慮到我們可以強行讓 \(x\) 取到某個值（比如 \(x_i\)）時，它所對應的值為 \(y_i\) 那麼可以將式子變成 \(f(x)=y_ik\) 那麼如何構造 \(k\) 便是我們要考慮的了，當 \(x\) 不取 \(x_i\)、在取其它 \(x_j\) 的時候，我們可以透過一個 \(x - x_j\) 來將其變為 \(0\)，可當 \(x_i\) 時，會有一些其他的東西混進去 (\(x-x_j\) 帶來的) 所以我們在構造一個 \(a_i-a_j\) 就可以將前面帶了的給去掉，考慮到上面是一個 \(0\) 就變成全 \(0\) 所以上面用乘法，所以將前面的去掉的部分我們也用乘法，兩個部分一除，便是答案（哇，好抽象）：

\(f(x)=\sum y_i\prod\frac{x-x_i}{x_i-x_j}\)

我們也可以將其理解為構造出來一個 \(k\) 來將函式扭成平滑曲線（可以從分段函式，到用平滑曲線來連結）

CF622F The Sum of the k-th Powers

首先我們要證明 \(\sum i^k\) 是個 \(k + 1\) 次多項式，然後這裡可以用 \(n^k-(n-1)^k\) 來證明差分後的次數會 \(-1\) 然後我們便可以證明上述命題（說了跟說了似的，算了自行推導），然後我們便可以考慮拉插，考慮到這些 \(x\) 值連續，所以我們便可以透過一些奇怪的做法將整個最佳化到 \(O(n)\)（哇，我真是個天才，這個也還不如自行推導）

#include <iostream>

using namespace std;
using ll = long long;

const int MaxN = 1e6 + 10, mod = 1e9 + 7;

ll f[MaxN], pr[MaxN], pl[MaxN], y[MaxN], n, k, ans, sum = 1;

ll qpow(ll a, ll b) {
  ll res = 1;
  for (ll i = 1; i <= b; i <<= 1) {
    (b & i) && (res = res * a % mod);
    a = a * a % mod;
  }
  return res;
}

int main() {
  cin >> n >> k, f[0] = pl[0] = pr[k + 3] = 1;
  for (int i = 1; i <= k + 2; i++) {
    f[i] = f[i - 1] * i % mod, sum = sum * (n - i) % mod;
    pl[i] = pl[i - 1] * (n - i) % mod;
  }
  for (int i = k + 2; i >= 1; i--) {
    pr[i] = pr[i + 1] * (n - i) % mod;
  }
  for (int i = 1; i <= k + 2; i++) {
    y[i] = (y[i - 1] + qpow(i, k)) % mod;
  }
  for (int i = 1; i <= k + 2; i++) {
    ll a = pl[i - 1] * pr[i + 1] % mod, b = f[i - 1] * f[k + 2 - i] % mod * ((k + 2 - i) % 2 ? -1 : 1) % mod;
    ans = (ans + a * qpow(b, mod - 2) % mod * y[i] % mod) % mod;
  }
  cout << (ans % mod + mod) % mod << endl;
  return 0;
}

P4463 [集訓隊互測 2012] calc

這個首先得出 dp 式子（什麼，這個不知道你學拉插幹什麼？）：

\(f_{i,j}=f{i-1,j-1}*j+f_{i,j - 1}\)

設 \(g_i=f_{n,i}\)，然後同樣用差分可以得到 \(f_{i,j} - f_{i,j - 1} = g_{i} - 1\)，然後看 dp 式子分析多項式次數：

\(g_i - 1 = g_{i - 1} + 1\)

那麼我們可以得到 \(g_i = g_{i - 1) + 2\)

所以 \(g_n = 2 * n\)

所以我們便知道 f_{n,i} 是一個 \(2 * n\) 次多項式，所以我們可以求出前面 \(2 * n + 1\) 個點的點值，然後拉插啟動！

#include <iostream>

using namespace std;
using ll = long long;

const int MaxN = 2e3 + 10;

ll f[MaxN][MaxN], a[MaxN], b[MaxN], n, k, p, mod, ans, sum = 1;

ll qpow(ll a, ll b) {
  ll res = 1;
  for (ll i = 1; i <= b; i <<= 1) {
    (b & i) && (res = res * a % mod);
    a = a * a % mod;
  }
  return res;
}

int main() {
  cin >> k >> n >> p, mod = p;
  for (int i = 0; i <= 2 * n + 1; i++) {
    f[0][i] = 1;
  }
  for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= 2 * n + 1; j++) {
      f[i][j] = (f[i - 1][j - 1] * j % mod + f[i][j - 1]) % mod;
    }
  }
  for (int i = 1; i <= 2 * n + 1; i++) {
    a[i] = i;
  }
  for (int i = 1; i <= 2 * n + 1; i++) {
    b[i] = f[n][i];
    sum = 1;
    for (int j = 1; j <= 2 * n + 1; j++) {
      if (i != j) sum = ((sum * (k - a[j]) % mod) + mod) % mod, sum = sum * (qpow(((a[i] - a[j]) % mod + mod) % mod, mod - 2)) % mod;
    }
    ans = (ans + sum * b[i] % mod) % mod;
  }
  for (int i = 1; i <= n; i++) {
    ans = ans * i % mod;
  }
  cout << ans << endl;
  return 0;
}

相關文章