拉格朗日插值与差分学习笔记（维护中）
概述
点值表达
对于多项式的表达有系数表达和点值表达。
点值表达：给定满足 \(y = f(x)\) 的 \(n\) 个点 \((x _ {1, \cdots, n}, y _ {1, \cdots, n})\) 可以表示出一个 \(n - 1\) 次多项式。
用途 - 1
多项式乘法。
目前没学到，后面填坑。
拉格朗日插值
本部分一切枚举下标都从 \(0\) 开始。
普通拉格朗日插值
插值是一种通过通过点值还原多项式的过程。拉格朗日插值是 OI 中较为常用的一种插值。
考虑构造 \(n\) 个函数，满足 $ f_i(x_i) = 1, \forall j \not = i, f_i(x_j) = 0 $。
有构造：
证明显然，对于 \(x = x_j, j \not = i\)，分子一定存在一个因式为 \(0\)。而对于 \(x = x_i\)，连乘的每一项都是 \(1\)。
那么则有：
暴力计算 \(f(x)\) 复杂度 \(\mathcal{O}(n ^ 2)\)，加上求逆元的普遍复杂度 \(\mathcal{O}(\log V)\)，总时间复杂度 \(\mathcal{O}(n ^ 2 \log V)\)。
重心拉格朗日插值
考虑优化。
考虑变形：
定义：
对每个 \(i\)，\(\mathcal{O}(n)\) 预处理 \(w_i\)，预处理总时间复杂度 \(\mathcal{O}(n^2)\)。
对每个 \(x\)，\(\mathcal{O}(n)\) 计算 \(\prod \limits _ {j} {(x - x_j)}\)，\(\mathcal{O}(n)\) 计算 \(\sum \limits _i \frac {w_i} {(x - x_i)}\)，单次计算 \(f(x)\) 的总时间复杂度 \(\mathcal{O}(n)\)。
代码实现
点击查看代码
#define int ll
const int N = 2e3 + 5;
const int mod = 998244353;
int qpow(int a, int b){
int res = 1;
while(b){
if(b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int inv(int x){
x = (x % mod + mod) % mod;
return qpow(x, mod - 2);
}
int n, k;
struct lagrange{
int n; // 多项式次数
struct node{
int x, y;
} p[N];
int w[N];
void init(int n_){
n = n_;
rep(i, 0, n){
w[i] = (p[i].y % mod + mod) % mod;
rep(j, 0, n){
if(i == j) continue;
w[i] = w[i] * inv(p[i].x - p[j].x) % mod;
}
}
}
int calc(int x){
int prod = 1, sum = 0;
rep(i, 0, n){
prod = prod * (x - p[i].x + mod) % mod;
sum = (sum + w[i] * inv(x - p[i].x) % mod + mod) % mod;
}
return prod * sum % mod;
}
} f;
void solve_test_case(){
n = read(), k = read();
rep(i, 0, n - 1){
f.p[i].x = read(), f.p[i].y = read();
}
f.init(n - 1);
write(f.calc(k));
}
自变量连续清况
实际题目中，常常出现点值自变量 \(x\) 连续的情况，针对这种情况有 \(\mathcal{O}(n)\) 复杂度的优化。
变为：
对分子，维护 \(x - j\) 的前缀积和后缀积，即可快速计算。即：
对分母，\(0 \le j \le n\)，实际上是这样的一个序列的全局积：\(i, i-1, \cdots, 2, 1, -1, -2, \cdots, -(n - i)\)。
拆成前后两段，分别为 \(i!\) 和 \((-1)^{n-i}(n-i)!\)。
最终结果为：
预处理阶乘，时间复杂度 \(\mathcal{O}(n)\)。
代码实现
点击查看代码
struct lagrange{
int n; // n 次多项式，[0, n] 为点值 x 坐标
int p[N], fac[N], inv[N], pre[N], suf[N];
void init(int n_){
n = n_;
fac[0] = 1;
rep(i, 1, n){
fac[i] = fac[i - 1] * i % mod;
}
inv[n] = qpow(fac[n], mod - 2);
per(i, n - 1, 0){
inv[i] = inv[i + 1] * (i + 1) % mod;
}
p[0] = 0;
}
int calc(int x){
if(x <= n) return p[x];
x %= mod;
pre[0] = x % mod;
rep(i, 1, n){
pre[i] = pre[i - 1] * (x - i + mod) % mod;
}
suf[n] = (x - n + mod) % mod;
per(i, n - 1, 0){
suf[i] = suf[i + 1] * (x - i + mod) % mod;
}
int res = 0;
rep(i, 0, n){
int tmp = p[i] * inv[n - i] % mod * inv[i] % mod;
if(i > 0) tmp = tmp * pre[i - 1] % mod;
if(i < n) tmp = tmp * suf[i + 1] % mod;
if((n - i) & 1) tmp = mod - tmp;
res = ((res + tmp) % mod + mod) % mod;
}
return res;
}
};
差分
定义多项式的差分，\(\Delta _ s ^ n f(x)\) 表示 \(f(x)\) 步长为 \(s\) 的 \(n\) 阶差分。
有：
用途 - 1
在实际问题中，往往需要自己分析所求多项式的次数，进而知道需要多少个点值来确定它。
结论：若多项式 \(f(x)\) 的 \(n\) 阶差分为常数，则 \(f(x)\) 的次数为 \(n\)，需要 \(n + 1\) 个点插值求出。
证明：
系数表达 \(f(x) = \sum _ {i = 0} a_i x ^ i\)。
那么 \(f(x + s)\) 与 \(f(x)\) 的 \(n\) 次项相减为：
\[a_n \left( (x + s) ^ n - x ^ n \right) \]其 \(n\) 次项系数显然为 \(0\)。
故每做一次差分，多项式次数都会 - 1。
例题
BZOJ 3453 | tyvj 1858 XLkxc
题意
给定 \(k, a, s, d, p = 1234567891\)。
求 \(h(s)\)。
思路
使用拉格朗日插值求解出 \(f, g, h\) 三个多项式。
首先，对于 \(f(x)\)，\(\Delta _ 1 ^ 1 f(x) = (x + 1) ^ k\) 为 \(k\) 次多项式，故 \(f(x)\) 为 \(k + 1\) 次多项式。
类似地，\(\Delta _ 1 ^ 1 g(x) = f(x + 1)\) 为 \(k + 1\) 次多项式，故 \(g(x)\) 为 \(k + 2\) 次多项式。
\(\Delta _ 1 ^ 1 h(x) = g(a + (x + 1)d)\) 为 \(k + 2\) 次多项式，故 \(h(x)\) 为 \(k + 3\) 次多项式。
代码实现中，我们求出每个函数自变量范围 \([0, n]\) 的函数值进行插值。应用自变量连续的拉格朗日插值，时间复杂度 \(\mathcal{O}(n)\)。
代码
点击查看代码
#define int ll
const int N = 130;
const int mod = 1234567891;
int qpow(int a, int b){
int res = 1;
while(b){
if(b & 1) res = res * a % mod;
a = a * a % mod;
b >>= 1;
}
return res;
}
int k, a, s, d;
struct lagrange{
int n; // n 次多项式，[0, n] 为点值 x 坐标
int p[N], fac[N], inv[N], pre[N], suf[N];
void init(int n_){
n = n_;
fac[0] = 1;
rep(i, 1, n){
fac[i] = fac[i - 1] * i % mod;
}
inv[n] = qpow(fac[n], mod - 2);
per(i, n - 1, 0){
inv[i] = inv[i + 1] * (i + 1) % mod;
}
p[0] = 0;
}
int calc(int x){
if(x <= n) return p[x];
x %= mod;
pre[0] = x % mod;
rep(i, 1, n){
pre[i] = pre[i - 1] * (x - i + mod) % mod;
}
suf[n] = (x - n + mod) % mod;
per(i, n - 1, 0){
suf[i] = suf[i + 1] * (x - i + mod) % mod;
}
int res = 0;
rep(i, 0, n){
int tmp = p[i] * inv[n - i] % mod * inv[i] % mod;
if(i > 0) tmp = tmp * pre[i - 1] % mod;
if(i < n) tmp = tmp * suf[i + 1] % mod;
if((n - i) & 1) tmp = mod - tmp;
res = ((res + tmp) % mod + mod) % mod;
}
return res;
}
} f, g, h;
void solve_test_case(){
k = read(), a = read(), s = read(), d = read();
f.init(k + 1);
f.p[0] = 0;
rep(i, 1, k + 1){
f.p[i] = (f.p[i - 1] + qpow(i, k)) % mod;
}
g.init(k + 2);
g.p[0] = 0;
rep(i, 1, k + 2){
g.p[i] = (g.p[i - 1] + f.calc(i)) % mod;
}
h.init(k + 3);
h.p[0] = g.calc(a);
rep(i, 1, k + 3){
h.p[i] = (h.p[i - 1] + g.calc(a + i * d)) % mod;
}
write(h.calc(s));
}
浙公网安备 33010602011771号