拉格朗日插值法

简介

众所周知，nnn 个点 (xi,yi)(x_i, y_i)(xi,yi)（任意两个点横坐标不相等）可以确定一个 n−1n - 1n−1 次多项式函数 y=f(x)y = f(x)y=f(x)。拉格朗日插值法可以根据这 nnn 个点求出这个多项式 f(x)f(x)f(x)。当然，实际应用中通常求出横坐标为 kkk 的点在该（ nnn 个点确定的）多项式函数上对应的纵坐标的值，代码实现中我们也只考虑这一问题。

一个直观的想法是利用待定系数法设 f(x)=an−1xn−1+⋯+a0x0f(x) = a_{n- 1} x ^ {n - 1} + \cdots + a_0 x ^ 0f(x)=an−1xn−1+⋯+a0x0，然后带入 nnn 个点得到一个 nnn 元一次方程组，然后用高斯消元解得系数。但这个方法和拉格朗日插值法相比有两个问题：一是时间复杂度为 O(n3)O(n^3)O(n3)，而拉格朗日法时间复杂度为 n2n^2n2；二就是系数可能解出小数，还可能很大，而拉格朗日法可以支持取模，并跳过“系数”这一中间步骤，直接求值。

拉格朗日插值法

我们以二次函数为例，看一看拉格朗日插值法的具体过程：已知 333 个点 (x1,y1)(x_1, y_1)(x1,y1)，(x2,y2)(x_2, y_2)(x2,y2)，(x3,y3)(x_3, y_3)(x3,y3)，求 f(x)f(x)f(x)。

拉格朗日（Joseph-Louis Lagrange，1736 ~ 1813）的做法非常巧妙地避开了解多元方程的过程：
令 f1(x)f_1(x)f1(x) 表示经过点 (x1,1)(x_1, 1)(x1,1)，(x2,0)(x_2, 0)(x2,0)，(x3,0)(x_3, 0)(x3,0) 的二次函数；
f2(x)f_2(x)f2(x) 表示经过点 (x1,0)(x_1, 0)(x1,0)，(x2,1)(x_2, 1)(x2,1)，(x3,0)(x_3, 0)(x3,0) 的二次函数；
f3(x)f_3(x)f3(x) 表示经过点 (x0,1)(x_0, 1)(x0,1)，(x2,0)(x_2, 0)(x2,0)，(x3,1)(x_3, 1)(x3,1) 的二次函数。
那么 f(x)=y1⋅f1(x)+y2⋅f2(x)+y3⋅f3(x)f(x) = y_1 \cdot f_1(x) + y_2 \cdot f_2(x) + y_3 \cdot f_3(x)f(x)=y1⋅f1(x)+y2⋅f2(x)+y3⋅f3(x)。

原因很简单，每个子函数确保经过一个点而不经过另外两个点。

而子函数的求法很简单，以 f1(x)f_1(x)f1(x) 为例：
f1(x)=0f_1(x) = 0f1(x)=0 的两根为 x=x2x = x_2x=x2 和 x=x3x = x_3x=x3，于是设 f1(x)=k(x−x2)(x−x3)f_1(x) = k (x - x_2) (x - x_3)f1(x)=k(x−x2)(x−x3)，再带入点 (x1,1)(x_1, 1)(x1,1)，得到 k=1(x1−x2)(x1−x3)k = \frac{1}{(x_1 - x_2)(x_1 - x_3)}k=(x1−x2)(x1−x3)1，于是 f1(x)=(x−x2)(x−x3)(x1−x2)(x1−x3)f_1(x) = \frac{(x - x_2) (x - x_3)}{(x_1 - x_2)(x_1 - x_3)}f1(x)=(x1−x2)(x1−x3)(x−x2)(x−x3)。

求高次函数与求二次函数的方法同理，可得
fi(x)=∏1≤j≤n,j≠i(x−xj)(xi−xj)f(x)=∑1≤i≤nfi(x)\begin{aligned} f_i(x) &= \prod_{1 \leq j \leq n, j \neq i} \frac{(x - x_j)}{(x_i - x_j)} \\ f(x) &= \sum_{1 \leq i \leq n} f_i(x) \end{aligned} fi(x)f(x)=1≤j≤n,j=i∏(xi−xj)(x−xj)=1≤i≤n∑fi(x)
于是，想求 f(k)f(k)f(k) 的值，将 kkk 代入上式即可，时间复杂度 O(n2)O(n^2)O(n2)（nnn 为次数）。

模板

洛谷 P4781 【模板】拉格朗日插值

#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 看作一个函数，即令fi(j)=dp[i][j]f_i(j) = dp[i][j]fi(j)=dp[i][j]（注意这个 fi(j)f_i(j)fi(j) 与上文中的“子函数”没有关系）那么，如果我们要求的 dp[i][j]dp[i][j]dp[i][j] 中的 jjj 值很大（例如 j=109j = 10^9j=109），我们就可以只计算 dp[i][1],dp[i][2],⋯,dp[i][p+1]dp[i][1], dp[i][2], \cdots, dp[i][p + 1]dp[i][1],dp[i][2],⋯,dp[i][p+1]（ppp 为 fi(x)f_i(x)fi(x) 的次数），并用点 (1,dp[i][1])(1, dp[i][1])(1,dp[i][1])，(2,dp[i][2])(2, dp[i][2])(2,dp[i][2])，…，(p+1,dp[i][p+1])(p + 1, dp[i][p + 1])(p+1,dp[i][p+1]) 确定多项式 fi(x)f_i(x)fi(x)，并快速求得 fi(j)f_i(j)fi(j)，即 dp[i][j]dp[i][j]dp[i][j]，时间复杂度为 O(p2)O(p^2)O(p2)。

这类优化的难点在于要准确地计算 ppp 的值，即 fi(x)f_i(x)fi(x) 的次数，接下来通过例题讲解如何计算 ppp。

例题一

洛谷 P4463 [集训队互测2012] calc

分析

发现我们只需要计算所有递增的合法序列的值之和，然后乘上 n!n!n! 即为答案，因为每种递增的合法序列任意打乱顺序仍然是合法的，并且原先就不同，打乱后也一定不同。

令 dp[i][j]dp[i][j]dp[i][j] 表示：长度为 iii 的所含元素值不超过 jjj 的递增的合法序列的值之和，考虑在第 iii 个位置放元素 jjj 还是放其他小于 jjj 的元素，本质即为一个背包问题，则dp[i][j]=j⋅dp[i−1][j−1]+dp[i][j−1]dp[i][j] = j \cdot dp[i - 1][j - 1] + dp[i][j - 1]dp[i][j]=j⋅dp[i−1][j−1]+dp[i][j−1]答案为 dp[n][k]dp[n][k]dp[n][k]，然后发现 k≤109k \leq 10^9k≤109，不可能直接 DP。

按照上文中的方法，我们令 fn(i)=dp[n][i]f_n(i) = dp[n][i]fn(i)=dp[n][i]，所求的就是 fn(k)f_n(k)fn(k)。接下来求出多项式 fn(x)f_n(x)fn(x) 的次数 ppp，然后我们就只需要 DP 出 dp[n][1]dp[n][1]dp[n][1] 到 dp[n][p+1]dp[n][p + 1]dp[n][p+1]，再用拉格朗日插值法就能算出 fn(k)f_n(k)fn(k) 了。

接下来推导 fn(x)f_n(x)fn(x) 的次数，令 g(n)g(n)g(n) 表示多项式 fn(x)f_n(x)fn(x) 的次数：
dp[i][j]=j⋅dp[i−1][j−1]+dp[i][j−1]fi(j)=j⋅fi−1(j−1)+fi(j−1)fi(j)−fi(j−1)=j⋅fi−1(j−1)\begin{aligned} dp[i][j] &= j \cdot dp[i - 1][j - 1] + dp[i][j - 1] \\ f_i(j) &= j \cdot f_{i - 1}(j - 1) + f_i(j - 1) \\ f_i(j) - f_i(j - 1) &= j \cdot f_{i - 1}(j - 1) \end{aligned} dp[i][j]fi(j)fi(j)−fi(j−1)=j⋅dp[i−1][j−1]+dp[i][j−1]=j⋅fi−1(j−1)+fi(j−1)=j⋅fi−1(j−1)设 fi(x)=∑i=0g(n)aixif_i(x) = \sum\limits_{i = 0}^{g(n)} a_i x ^ifi(x)=i=0∑g(n)aixi，将 jjj 和 j−1j - 1j−1 暴力代入 fi(j)−fi(j−1)f_i(j) - f_i(j - 1)fi(j)−fi(j−1) 这个式子，发现 ag(i)jg(n)a_{g(i)} j^{g(n)}ag(i)jg(n) 这个最高次项被消掉了（代入后有关最高次项的部分仅为 ag(i)jg(i)−ag(i)(j−1)g(i)a_{g(i)} j^{g(i)} - a_{g(i)} (j - 1)^{g(i)}ag(i)jg(i)−ag(i)(j−1)g(i)）！

于是得到 fi(j)−fi(j−1)f_i(j) - f_i(j - 1)fi(j)−fi(j−1) 的次数为 g(i)−1g(i) - 1g(i)−1，又因为 j⋅fi−1(j−1)j \cdot f_{i - 1}(j - 1)j⋅fi−1(j−1) 的次数为 g(i−1)+1g(i - 1) + 1g(i−1)+1，所以
g(i)−1=g(i−1)+1g(i)=g(i−1)+2\begin{aligned} g(i) - 1 &= g(i - 1) + 1\\ g(i) &= g(i - 1) + 2 \end{aligned} g(i)−1g(i)=g(i−1)+1=g(i−1)+2又因为 g(0)=0g(0) = 0g(0)=0（f0(x)=dp[0][x]=1f_0(x) = dp[0][x] = 1f0(x)=dp[0][x]=1）所以 g(n)=2ng(n) = 2ng(n)=2n，证得 fn(x)f_n(x)fn(x) 的次数为 2n2n2n。

然后我们只需要用朴素的 DP 求得 dp[n][1]dp[n][1]dp[n][1]，dp[n][2]dp[n][2]dp[n][2]，…，dp[n][2n+1]dp[n][2n + 1]dp[n][2n+1]（注意点数要求比次数多一才能得到正确的多项式），并用拉格朗日插值法求得 dp[n][k]dp[n][k]dp[n][k] 即可。

代码

#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;
}

例题二

CF995F Cowmpany Cowmpensation

题意：给定整数 nnn 和 DDD（1≤n≤30001 \leq n \leq 30001≤n≤3000，1≤D≤1091 \leq D \leq 10^91≤D≤109）以及一个 nnn 个结点的树，要求给每个结点分配一个 [1,D][1, D][1,D] 之间的整数作为权值，并且满足父亲结点权值大于等于儿子结点，求方案总数。

分析

令 dp[u][i]dp[u][i]dp[u][i] 表示：以 uuu 为根的子树中，每个结点的权值都在 [1,i][1,i][1,i] 内的方案数，同样是一个背包
dp[u][i]=dp[u][i−1]+∑vis a son of udp[v][i−1]dp[u][i] = dp[u][i - 1] + \sum_{v \text{ is a son of } u} dp[v][i - 1]dp[u][i]=dp[u][i−1]+v is a son of u∑dp[v][i−1]g(n)g(n)g(n) 定义与上题类似，然后得到
g(u)−1=∑vis a son of ug(v)g(u)=∑vis a son of ug(v)+1\begin{aligned} g(u) - 1 &= \sum_{v \text{ is a son of } u} g(v)\\ g(u) &= \sum_{v \text{ is a son of } u} g(v) + 1 \end{aligned} g(u)−1g(u)=v is a son of u∑g(v)=v is a son of u∑g(v)+1注意边界 g(v)=[vis a leaf ]g(v) = [v \text{ is a leaf }]g(v)=[v is a leaf ]，因为对于一个叶子 uuu 有 dp[u][i]=idp[u][i] = idp[u][i]=i。因此这就是一个子树大小的 DP 式，于是 g(1)=ng(1) = ng(1)=n，暴力算得 dp[1][1]dp[1][1]dp[1][1]，dp[1][2]dp[1][2]dp[1][2]，…，dp[1][n+1]dp[1][n + 1]dp[1][n+1]，再拉格朗日即可。

代码

#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

上面两题的“点”的横坐标有个规律：是连续的 p+1p + 1p+1 个正整数。结合拉格朗日插值法的分子分母的特征，发现可以用前缀积和后缀积优化拉格朗日插值法的内层循环代码，使时间复杂度由 O(p2)O(p^2)O(p2) 优化为 O(p)O(p)O(p)，但是复杂度的瓶颈在于开头的朴素 DP，所以没有提这个方法。

C++ 拉格朗日插值法优化 DP相关推荐

洛谷 P4463 [集训队互测 2012] calc（拉格朗日插值优化DP）
整理的算法模板合集: ACM模板点我看算法全家桶系列!!! 实际上是一个全新的精炼模板整合计划 Weblink https://www.luogu.com.cn/problem/P4463 Prob ...
【BZOJ2655】calc，dp+拉格朗日插值法
传送门思路: 好题这个题目做法非常多,我知道的起码还有倍增法和预处理伯努利数这里仅说一下我的想法从原始dp的思路出发 f[i][j]表示i个元素中选取j个的分数写出dp方程 f[i][j]= ...
rstudio拉格朗日插值法_电力窃漏电用户识别案例
一.案例综述案例编号:102003 案例名称:电力.热力.燃气及水生产和供应业--电力窃漏电用户识别作者姓名(或单位.或来源):朱江案例所属行业:D442 电力供应案例所用软件:R 案例包含知 ...
拉格朗日插值法及应用
拉格朗日插值法一般方法重心拉格朗日插值法应用 bzoj4559:成绩比较 bzoj2655: calc bzoj3453:XLkxc 拉格朗日插值法快速根据点值逼近函数在取点大于nnn的情况 ...
rstudio拉格朗日插值法_拉格朗日插值法学习笔记
拉格朗日插值法是一个根据点对求回原函数的算法,原理挺好懂的. 原理和优化方法上面的大佬都讲得很好. 其实主要就是这个式子: 然后暴力算这个式子的话是每求一项f(k)的时间复杂度都是n^2. 这个时间很 ...
拉格朗日插值法求多项式系数（附代码）
写在前面: 学了拉格朗日插值法之后发现大家都说可以在O(n^2)时间内得到多项式系数,但是没有找到代码,网上找了很多资料又因为我太弱了没能看懂,最后在emofunx学长的帮助下终于搞明白了. 由于太弱 ...
拉格朗日插值法学习笔记
文章目录 STEP1 问题引入 STEP2 拉格朗日插值公式存在性唯一性代码模板 STEP3 算法优化重心拉格朗日插值 xxx连续取值的情况 STEP4 经典例题 The Sum of the ...
rstudio拉格朗日插值法_拉格朗日插值法
我们遇到的问题是,给定一个多项式的点值表示和一个数,求出这个数带入多项式后的值. 这个问题如果用待定系数法,可以使用高斯消元,但是复杂度是$O(n^3)$的,无法通过本题. 所以我们来引入拉格朗日 ...
[学习笔记]矩阵乘法及其优化dp
1.定义: $c[i][j]=\sum a[i][k]\times b[k][j]$ 所以矩阵乘法有条件,(n*m)*(m*p)=n*p 即第一个矩阵的列数等于第二个矩阵的行数,否则没有意义. 2.结 ...

C++ 拉格朗日插值法优化 DP

文章目录

拉格朗日插值法

简介