拉格朗日插值法
简介
众所周知, 个点 (任意两个点横坐标不相等)可以确定一个 次多项式函数 。拉格朗日插值法可以根据这 个点求出这个多项式 。当然,实际应用中通常求出横坐标为 的点在该( 个点确定的)多项式函数上对应的纵坐标的值,代码实现中我们也只考虑这一问题。
一个直观的想法是利用待定系数法设 ,然后带入 个点得到一个 元一次方程组,然后用 高斯消元 解得系数。但这个方法和拉格朗日插值法相比有两个问题:一是时间复杂度为 ,而拉格朗日法时间复杂度为 ;二就是系数可能解出小数,还可能很大,而拉格朗日法可以支持取模,并跳过“系数”这一中间步骤,直接求值。
拉格朗日插值法
我们以二次函数为例,看一看拉格朗日插值法的具体过程:已知 个点 , , ,求 。
拉格朗日(Joseph-Louis Lagrange,1736 ~ 1813)的做法非常巧妙地避开了解多元方程的过程:
令
表示经过点
,
,
的二次函数;
表示经过点
,
,
的二次函数;
表示经过点
,
,
的二次函数。
那么
。
原因很简单,每个子函数确保经过一个点而不经过另外两个点。
而子函数的求法很简单,以
为例:
的两根为
和
,于是设
,再带入点
,得到
,于是
。
求高次函数与求二次函数的方法同理,可得
于是,想求
的值,将
代入上式即可,时间复杂度
(
为次数)。
模板
#include <bits/stdc++.h>
const int MOD = 998244353;
const int MAXN = 2000;
int Mul(const int &a, const int &b) {
return (long long)a * b % MOD;
}
int Inv(int x) {
int y = MOD - 2, ret = 1;
while (y) {
if (y & 1)
ret = Mul(ret, x);
x = Mul(x, x);
y >>= 1;
}
return ret;
}
int N, K, X[MAXN + 5], Y[MAXN + 5];
int main() {
scanf("%d%d", &N, &K); // 求 f(K)
for (int i = 1; i <= N; i++)
scanf("%d%d", &X[i], &Y[i]);
int Ans = 0;
for (int i = 1; i <= N; i++) {
int x = Y[i], y = 1;
for (int j = 1; j <= N; j++)
if (j != i) {
x = Mul(x, (K - X[j] + MOD) % MOD);
y = Mul(y, (X[i] - X[j] + MOD) % MOD);
}
Ans = (Ans + Mul(x, Inv(y))) % MOD;
}
printf("%d", Ans);
}
DP 优化
思路
如果没有接触过可能很难想到这个与 DP 的联系。事实上,我们可以将某一维的 DP 看作一个函数,即令 (注意这个 与上文中的“子函数”没有关系)那么,如果我们要求的 中的 值很大(例如 ),我们就可以只计算 ( 为 的次数),并用点 , ,…, 确定多项式 ,并快速求得 ,即 ,时间复杂度为 。
这类优化的难点在于要准确地计算 的值,即 的次数,接下来通过例题讲解如何计算 。
例题一
分析
发现我们只需要计算所有递增的合法序列的值之和,然后乘上 即为答案,因为每种递增的合法序列任意打乱顺序仍然是合法的,并且原先就不同,打乱后也一定不同。
令 表示:长度为 的所含元素值不超过 的递增的合法序列的值之和,考虑在第 个位置放元素 还是放其他小于 的元素,本质即为一个背包问题,则 答案为 ,然后发现 ,不可能直接 DP。
按照上文中的方法,我们令 ,所求的就是 。接下来求出多项式 的次数 ,然后我们就只需要 DP 出 到 ,再用拉格朗日插值法就能算出 了。
接下来推导
的次数,令
表示多项式
的次数:
设
,将
和
暴力代入
这个式子,发现
这个最高次项被消掉了(代入后有关最高次项的部分仅为
)!
于是得到
的次数为
,又因为
的次数为
,所以
又因为
(
)所以
,证得
的次数为
。
然后我们只需要用朴素的 DP 求得 , ,…, (注意点数要求比次数多一才能得到正确的多项式),并用拉格朗日插值法求得 即可。
代码
#include <bits/stdc++.h>
const int MAXN = 500;
int N, K, P;
int Dp[MAXN + 5][2 * MAXN + 1 + 5];
int Add(int a, const int &b) {
a += b; return (a >= P) ? (a - P) : a;
}
int Mul(const int &a, const int &b) {
return (long long)a * b % P;
}
int Inv(int x) {
int y = P - 2, ret = 1;
while (y) {
if (y & 1)
ret = Mul(ret, x);
x = Mul(x, x);
y >>= 1;
}
return ret;
}
int main() {
scanf("%d%d%d", &K, &N, &P);
int M = 2 * N + 1;
for (int i = 0; i <= M; i++)
Dp[0][i] = 1;
for (int i = 1; i <= N; i++)
for (int j = i; j <= M; j++)
Dp[i][j] = Add(Dp[i][j - 1], Mul(Dp[i - 1][j - 1], j));
int Ans = 0, Fac = 1;
for (int i = 1; i <= N; i++)
Fac = Mul(Fac, i);
for (int i = 1; i <= M; i++) {
int x = Dp[N][i], y = 1;
for (int j = 1; j <= M; j++)
if (i != j) {
x = Mul(x, (K >= j) ? (K - j) : (K - j + P));
y = Mul(y, (i >= j) ? (i - j) : (i - j + P));
}
Ans = Add(Ans, Mul(x, Inv(y)));
}
printf("%d", Mul(Ans, Fac));
return 0;
}
例题二
题意:给定整数 和 ( , )以及一个 个结点的树,要求给每个结点分配一个 之间的整数作为权值,并且满足父亲结点权值大于等于儿子结点,求方案总数。
分析
令
表示:以
为根的子树中,每个结点的权值都在
内的方案数,同样是一个背包
定义与上题类似,然后得到
注意边界
,因为对于一个叶子
有
。因此这就是一个子树大小的 DP 式,于是
,暴力算得
,
,…,
,再拉格朗日即可。
代码
#include <bits/stdc++.h>
const int MAXN = 3000;
const int MOD = 1000000007;
int N, D, M;
std::vector<int> G[MAXN + 5];
int Dp[MAXN + 5][MAXN + 5];
int Add(int a, const int &b) {
a += b; return (a >= MOD) ? (a - MOD) : a;
}
int Mul(const int &a, const int &b) {
return (long long)a * b % MOD;
}
int Inv(int x) {
int y = MOD - 2, ret = 1;
while (y) {
if (y & 1)
ret = Mul(ret, x);
x = Mul(x, x);
y >>= 1;
}
return ret;
}
void Dfs(int u) {
for (int v: G[u])
Dfs(v);
for (int i = 1; i <= M; i++) {
int tmp = 1;
for (int v: G[u])
tmp = Mul(tmp, Dp[v][i]);
Dp[u][i] = Add(Dp[u][i - 1], tmp);
}
}
int main() {
scanf("%d%d", &N, &D);
for (int i = 2; i <= N; i++) {
int p; scanf("%d", &p);
G[p].push_back(i);
}
M = N + 1;
Dfs(1);
int Ans = 0;
for (int i = 1; i <= M; i++) {
int x = Dp[1][i], y = 1;
for (int j = 1; j <= M; j++)
if (i != j) {
x = Mul(x, (D >= j) ? (D - j) : (D - j + MOD));
y = Mul(y, (i >= j) ? (i - j) : (i - j + MOD));
}
Ans = Add(Ans, Mul(x, Inv(y)));
}
printf("%d", Ans);
return 0;
}
A trick
上面两题的“点”的横坐标有个规律:是连续的 个正整数。结合拉格朗日插值法的分子分母的特征,发现可以用前缀积和后缀积优化拉格朗日插值法的内层循环代码,使时间复杂度由 优化为 ,但是复杂度的瓶颈在于开头的朴素 DP,所以没有提这个方法。