「学习笔记」拉格朗日插值

感谢 clf 指出笔误！/qq

一般形式的构造过程

\((n+1)\) 个横坐标不同的点可以唯一确定一个 \(n\) 次多项式，假如知道答案是个关于某个值的多项式以及其项数和若干个点，那么就可以拉格朗日插值出多项式，算出答案。

简单来讲，我们需要构造一个 \(n\) 次多项式，使得它满足在已知的点 \(x_i\) 处取值为 \(y_i\)。

构造 \(f_i(x)\)，使得其仅在 \(x_i\) 处取值为 \(1\)，其他 \(x_j,j\neq i\) 处取值为 \(0\)。

那么 \(y_if_i(x)\) 仅在 \(x_i\) 处取值为 \(y_i\)，其他 \(x_j,j\neq i\) 处取值为 \(0\)。

对 \(y_if_i(x)\) 求和即得到了我们要求的多项式：在 \(x_i\) 处取值为 \(y_i\) 的 \(n\) 次多项式。

如何构造 \(f_i(x)\)？

• 如果要在 \(x_j,j\neq i\) 处取值为 \(0\)，可以让 \((x-x_j),j\neq i\) 乘起来，即为 \(\prod_{j\neq i}(x-x_j)\)，为描述方便记作 \(g(x)\)，基于这个式子，不管怎么对其乘除，在 \(x_j,j\neq i\) 处取值都为 \(0\)。

• 如果要在 \(x_i\) 处取值为 \(1\)，则让 \(f_i(x)=\frac{g(x)}{g(x_i)}\)，这样当 \(x=x_i\) 时，\(f_i(x_i)=\frac{g_{x_i}}{g_{x_i}}=1\)，注意到 \(x_i\neq j,i\neq j\)，所以不会出现除 \(0\) 的情况。

综上所述，我们构造出了 \(f_i(x)=\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}\)，题目中所求的 \(n\) 次多项式即为：

\[\begin{aligned} &\sum _{i=1}^{n}y_if_i(x) \\ &\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \end{aligned} \]

若 \(x_i\) 是连续的 \(n\) 个整数的 \(\mathcal{O}(n)\) 做法

代入值到 \(x\) 中直接对这个式子计算，每一次枚举 \(i\) 后，最后一块算分母的逆元，只会算 \(n\) 次逆元，时间复杂度为 \(\mathcal{O}(n^2)\)。

特别地，若 \(x_i\) 是连续的 \(n\) 个整数，则可以 \(\mathcal{O}(n)\) 计算。

设 \(pre_i=\prod_{j=1}^i x-j,suf_i=\prod_{j=i+1}^nx-j\)，也就是分子的前缀/后缀积，\(fac_i=i!\)。

则

\[\begin{aligned} f(x)&=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \\ &=\sum_{i=1}^ny_i\frac{pre_{i-1}\cdot suf_{i+1}}{fac_{i-1}\cdot fac_{n-i}\cdot (-1)^{n-i}} \end{aligned} \]

可以线性递推求逆元，\(\mathcal{O}(n+\log mod)\) 求逆元也是可以接受的，故复杂度为 \(\mathcal{O}(n)\)

重心拉格朗日插值法

如果每加一次点就要询问，如何解决？推一下式子试试。

\[\begin{aligned} f(x)&=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \\ &=\sum_{i=1}^ny_i\frac{\prod_{j\neq i}(x-x_j)}{\prod_{j\neq i}(x_i-x_j)} \\ &=\sum_{i=1}^ny_i\frac{\prod_{j=1}^n(x-x_j)}{\left(\prod_{j\neq i}(x_i-x_j)\right )\cdot (x-x_i)} \\ &=\prod_{j=1}^n(x-x_j)\sum_{i=1}^n\frac{y_i}{\left(\prod_{j\neq i}(x_i-x_j)\right )\cdot (x-x_i)} \end{aligned} \]

设 \(g=\prod_{j=1}^n(x-x_j),w(i)=\prod_{j\neq i}^n(x_i-x_j)\)

则原式 \(=g\sum_{i=1}^n\frac{y_i}{w(i)\cdot (x-x_i)}\).

每加入一个点，\(\mathcal{O}(n)\) 计算出 \(g\) 并更新一下 \(w(i)\) 再 \(\mathcal{O}(n)\) 计算 \(f(x)\) 即可。

栗题一

ABC208F

给定 \(n,m,k\)，计算 \(f(n,m)\) 的值，模 \(10^9+7\)．

\[\begin{aligned} \displaystyle f(n, m)& = \begin{cases} 0 & (n = 0) \newline \\ n^K & (n \gt 0, m = 0) \newline \\ f(n-1, m) + f(n, m-1) & (n \gt 0, m \gt 0) \end{cases} \end{aligned} \]

\(0 \leq N \leq 10^{18},0 \leq M \leq 30,1 \leq K \leq 2.5 \times 10^6\)．

易发现答案为 \(1^k,2^k,3^k,\cdots ,n^k\) 的 \(m\) 阶前缀和的第 \(n\) 项值。

考虑 \(i^k\) 对答案的贡献次数，是对 \(f_0=1,f_k=0(k>0)\) 的数组作 \(m\) 阶前缀和后的第 \((n-i)\) 项，也就是 \(f_{n-i}\)。其值为 \(\binom{n-i+m-1}{m-1}\)．

从组合意义上考虑，作一阶前缀和为 \(1,1,1,\cdots\)，设后面第 \(i\) 阶前缀和的 \(f_j=g_{i,j}\)，有递推式 \(g_{i,j}=g_{i-1,j}+g_{i,j-1}\)，恰为只能向右、下走的格点计数。

故贡献的次数为 \((1,0)\to (m,n-i)\) 的路径数，即为 \(\binom{n-i+m-1}{m-1}\)。

所以答案为：

\[\sum_{i=0}^n i^k\binom{n-i+m-1}{m-1} \]

考虑 \(i^k\binom{n-i+m-1}{m-1}\) 中 \(i^k\) 为关于 \(i\) 的 \(k\) 次单项式。

\(\binom{n-i+m-1}{m-1}=\frac{(n-i+m-1)!}{(m-1)!(n-i)!}\)，其中 \((n-i+m-1)!\) 为关于 \(i\) 的 \((n-i+m-1)\) 次多项式，\((n-i)!\) 为关于 \(i\) 的 \((n-i)\) 次多项式，\((m-1)!\) 是常数。

故 \(i^k\binom{n-i+m-1}{m-1}\) 是关于 \(i\) 的 \((m+k-1)\) 次多项式。

对于 \(i\in [0,n]\)，对这个多项式求和即为答案。

其可以写成 \(ans=a_0s_0+a_1s_1+a_2s_2+\cdots+a_{m+k-1}s_{m+k-1}\) 的形式，其中 \(s_i=\sum\limits_{j=0}^ij^k\)．

由于 \(s_i\) 是关于 \(i\) 的 \((i+1)\) 次多项式，故答案为关于 \(n\) 的 \((m+k)\) 次多项式。

设 \(l=m+k+1\)，选出前 \(l\) 个连续的整数，算出其作 \(m\) 阶前缀和的答案，拉格朗日插值即可。

时间复杂度 \(\mathcal{O}(mk)\)．

#include<iostream>
#include<cstdio>
#include<algorithm>
typedef long long ll;
const ll mod = 1000000007;
template <typename T> T Max(T x, T y) { return x > y ? x : y; }
template <typename T> T Min(T x, T y) { return x < y ? x : y; }
template <typename T> T Add(T x, T y) { return (x + y >= mod) ? (x + y - mod) : (x + y); }
template <typename T> T Mod(T x) { return (x >= mod) ? (x - mod) : (x < 0 ? (x + mod) : x); }
template <typename T> T Mul(T x, T y) { return x * y % mod; }
template <typename T>
T &read(T &r) {
	r = 0; bool w = 0; char ch = getchar();
	while(ch < '0' || ch > '9') w = ch == '-' ? 1 : 0, ch = getchar();
	while(ch >= '0' && ch <= '9') r = (r << 3) + (r <<1) + (ch ^ 48), ch = getchar();
	return r = w ? -r : r;
}
ll qpow(ll x, ll y) { ll sumq = 1; while(y) { if(y & 1) sumq = sumq * x % mod; x = x * x % mod; y >>= 1; } return sumq; }
const int N = 2600100;
ll n;
int m, k;
ll a[N], fac[N], inv[N], pre[N], suf[N], ans;
signed main() {
	read(n); read(m); read(k); int l = m+k+1;
	if(!n) { puts("0");	return 0; }
	for(int i = 1; i <= l; ++i) a[i] = qpow(i, k);
	for(int j = 1; j <= m; ++j)
		for(int i = 1; i <= l; ++i)
			a[i] = Add(a[i], a[i-1]);
	inv[0] = fac[0] = 1; for(int i = 1; i <= l; ++i) fac[i] = fac[i-1] * i % mod;
	inv[l] = qpow(fac[l], mod-2); n %= mod;
	for(int i = l-1; i; --i) inv[i] = inv[i+1] * (i+1) % mod;
	pre[0] = 1; for(int i = 1; i <= l; ++i) pre[i] = pre[i-1] * (n - i) % mod;
	suf[l+1] = 1; for(int i = l; i; --i) suf[i] = suf[i+1] * (n - i) % mod;
	for(int i = 1; i <= l; ++i)
		ans = Add(ans, Mod(a[i] * pre[i-1] % mod * suf[i+1] % mod * inv[i-1] % mod * inv[l-i] % mod * (((l-i)&1) ? -1 : 1)));
	printf("%lld\n", ans);
	return 0;
}