拉格朗日插值法 (Lagrange interpolation approach) 学习笔记
Lagrange interpolation approach 是要解决一种如下的问题：
给定 \(n\) 个坐标，\((x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\)，确定一个多项式 \(f(x) = a_0 + a_1x + a_2x^2 + \dots + a_dx^d\) 满足：
\[f(x_1) = y_1
\]
\[f(x_2) = y_2
\]
\[\dots
\]
\[f(x_n) = y_n
\]
一、高斯消元
回想一下学 FFT 的时候，使用点值表示法，用 \(k + 1\) 个点表示一个 \(k\) 次的多项式。
那么我们可以设：
\[f(x) = a_0 + a_1x + a_2x^2 + \dots + a_{n - 1}x^{n - 1}
\]
然后我们可以列出方程组：
\[\begin{cases}a_0 + a_1x_1 + a_2x_1^2 + \dots + a_{n - 1}x_1^{n - 1} = y_1\\a_0 + a_1x_2 + a_2x_2^2 + \dots + a_{n - 1}x_2^{n - 1} = y_2\\\dots\\a_0 + a_1x_n + a_2x_n^2 + \dots + a_{n - 1}x_n^{n - 1} = y_n \end{cases}
\]
然后就是一次简单的高斯消元，时间复杂度 \(O(n ^ 3)\)。
二、拉格朗日插值
如果对推导感兴趣可以看这个：https://oi-wiki.org/math/poly/lagrange/
直接放公式：
\[f(x) = \sum_{i = 1}^{n} y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}
\]
三、洛谷 P4781 【模板】拉格朗日插值
https://www.luogu.com.cn/problem/P4781
#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/hash_policy.hpp>
using namespace std;
using namespace __gnu_pbds;
#ifdef LOCAL
#include "algo/debug.h"
#else
#define debug(...) 42
#endif
typedef long long ll;
typedef pair < int, int > PII;
typedef int itn;
mt19937 RND_MAKER (chrono :: steady_clock :: now ().time_since_epoch ().count ());
inline ll randomly (const ll l, const ll r) {return (RND_MAKER () ^ (1ull << 63)) % (r - l + 1) + l;}
#define int long long
const double pi = acos (-1);
//__gnu_pbds :: tree < Key, Mapped, Cmp_Fn = std :: less < Key >, Tag = rb_tree_tag, Node_Upadte = null_tree_node_update, Allocator = std :: allocator < char > > ;
//__gnu_pbds :: tree < PPS, __gnu_pbds :: null_type, less < PPS >, __gnu_pbds :: rb_tree_tag, __gnu_pbds :: tree_order_statistics_node_update > tr;
inline int read () {
int x = 0, f = 0;
char c = getchar ();
for ( ; c < '0' || c > '9' ; c = getchar ()) f |= (c == '-');
for ( ; c >= '0' && c <= '9' ; c = getchar ()) x = (x << 1) + (x << 3) + (c & 15);
return !f ? x : -x;
}
const int mod = 998244353;
const int N = 2e3 + 5;
int x[N], y[N], n, k, ans;
inline int pow_mod (int a, int b, int p) {
int res = 1;
while (b) {
if (b & 1) res = res * a % p;
b >>= 1;
a = a * a % p;
}
return res;
}
signed main () {
n = read (), k = read ();
for (int i = 1;i <= n; ++ i) {
x[i] = read (), y[i] = read ();
}
for (int i = 1;i <= n; ++ i) {
int tmp = y[i];
for (int j = 1;j <= n; ++ j) {
if (i != j) {
int l = k - x[j];
l = (l % mod + mod) % mod;
tmp = tmp * l % mod;
l = x[i] - x[j];
l = (l % mod + mod) % mod;
tmp = tmp * pow_mod (l, mod - 2, mod) % mod;
}
}
ans = (ans + tmp) % mod;
}
printf ("%lld\n", (ans % mod + mod) % mod);
return 0;
}
四、Codeforces 622F - The Sum of the k-th Powers
声明：博主数学很菜，所以并不会证明答案就是一个多项式 /kk /kk /kk。
我们把前 \(k + 2\) 个点都拎出来，比如第 \(i\) 个点就是 \(\displaystyle \left(i, \sum_{j = 1}^{i} j ^ k\right)\)。
然后就做一遍拉格朗日插值，把 \(f(n)\) 输出即可。
但是这样直接做是 \(O(k^2)\) 的，\(k\) 是 1e6 级别的，怎么办？
夹带力度！爆推式子：
我们暂且将 \(w_u\) 设为 \(\displaystyle \sum_{i = 1}^{u} i^k\)。
\[f(x) = \sum_{i = 1}^{n} w_i \prod_{j \neq i} \frac{x - j}{i - j}
\]
这个 \(j \neq i\) 的限制条件十分的烦人，我们可以拆成两个区间：
\[f(x) = \sum_{i = 1}^{n} \left(w_i \times \prod_{j = 1}^{i - 1} \frac{x - j}{i - j} \times \prod_{j = i + 1}^{n}\frac{x - j}{i - j}\right)
\]
\[f(x) = \sum_{i = 1}^{n} \left(w_i \times \frac{\prod\limits_{j = 1}^{i - 1} (x - j)}{(i - 1)!} \times \frac{\prod\limits_{j = i + 1}^{n} (x - j)}{(-1)^{n - i} \times (n - i)!}\right)
\]
特别的，当 \(k \leq 0\) 时，\(k! = 1\)。
分子可以直接前缀积算，然后分母可以预处理阶乘的逆元，就做完了。