拉格朗日插值学习笔记

\(\large{对于一个关于x的n次多项式f(x)，若知道其中的n+1个点，拉格朗日插值可以在\text{O}(\text{n}^2)的时间复杂度内求出f(k)。}\)
\(\\\)
\(\large{式子是这样的：\\ f\left( k\right) =\sum \limits^{n}_{i=0}y_{i_{j}}\prod \limits^{n}_{j=0,j\neq i}\dfrac {k-x_{j}}{x_{i}-x_{j}}\\关于正确性的证明，我还没理解。贴一下大佬的说法。}\)
\(\\\)

\(\\\)
\(\Large\textbf{Code: }\)

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

const int N = 2e3 + 5;
const int p = 998244353;

int n, k, x[N], y[N];

void read(int &x) {
	int flg = 1; x = 0;
	char c = getchar();
	while (!isdigit(c)) { if (c == '-') flg = -1; c = getchar(); }
	while (isdigit(c)) x = x * 10 + c - '0', c = getchar();
	x *= flg;
}

int pow(int x, int y) {
	ll a = x, ret = 1;
	for (; y ; y >>= 1, a = a * a % p) if (y & 1) ret = ret * a % p;
	return  ret;
}

int main() {
	read(n), read(k);
	for (int i = 1; i <= n; ++i) read(x[i]), read(y[i]);
	int s, f, ans = 0;
	for (int i = 1; i <= n; ++i) {
		f = y[i], s = 1;
		for (int j = 1; j <= n; ++j) {
			if (i == j) continue;
			f = 1ll * f * ((k - x[j] + p) % p) % p, s = 1ll * s * ((x[i] - x[j] + p) % p) % p;
		}
		s = pow(s, p - 2);
		ans = (ans + 1ll * s * f % p) % p;
	}
	printf("%d\n", ans);
	return 0;
} 

\(\\\)
\(\large{再说一下拓展，如何运用拉格朗提插值在\text{O}(\text{klogk})的时间内求出\sum \limits^{n}_{i=1}i^{k}，其中 n \leq 1e9, k \leq 1e6。}\)
\(\\\)
\(\large{首先证明一下\sum \limits^{n}_{i=1}i^{k}这个式子是一个k+1次的多项式。\\对于一个数列a \{n\}，把a \{n\}的相邻元素两两做差，得到b \{n\}，一般把数列b \{n\}成为a \{n\}的阶差数列。记数列b \{n\}为数列a \{n\}的一阶差分。\\如果再对数列b \{n\}进行差分，得到的就是数列a \{n\}的二阶差分，以此类推。\\定理1：数列 {a_n} 是一个 \text {p} 阶等差数列的充要条件是数列的通项 a_n 为 \text {n} 的一个\text {p} 次多项式。\\证明：\\有了上面那个定理，我们只关心最高项即可。\\设数列 { a_n} 的通项 a_n是一个关于 \text {n} 的 \text {v} 次多项式，即f\left( x\right) =\sum \limits^{\nu }_{i=0}u_{i}\cdot x^{i}。\\令数列b \{n\}为a \{n\}的一阶差分，那么b \{n\}的通项公式为f\left( x+1\right) -f\left( x\right) =\sum \limits^{v}_{i=0}u_{i}\times \left( x+1\right) ^{i}-\sum \limits^{v}_{i=0}u_{i}\times x_{i}。\\因为只考虑最高项，所以式子为u_v \times (x+1)^v - u_v \times x^v，然后用二项式定理把(x+1)^v展开，只考虑最高项，显然系数为1，那么式子就变成了u_v \times x^v - u_v \times x^v = 0。\\我们发现数列b \{n\}的通项公式中x^v被消掉了，所以数列b\{n\}是关于n的一个v-1次多项式。也就是说，做一次差分之后数列的通项公式的多项式次数会-1。\\那么数列a\{n\}在做p-1次差分后的到一个0次多项式，即数列a\{n\}的通项a_n为一个关于n的p次多项式。\\我们再考虑题目中的数列a\{n\}，即为\\ \sum \limits^{0}_{i=1}i^{k}, \sum \limits^{1}_{i=1}i^{k}, \sum\limits^{2}_{i=1}i^{k}... \sum\limits^{n}_{i=1}i^{k}，然后做差得到\\ \sum \limits^{1}_{i=1}i^{k}-\sum \limits^{0}_{i=1}i^{k}, \sum \limits^{2}_{i=1}i^{k}-\sum \limits^{1}_{i=1}i^{k}...\sum \limits^{n}_{i=1}i^{k}-\sum \limits^{n-1}_{i=1}i^{k}\\ 即1^k, 2^k, 3^k,...,n^k，显然这个数列的通项公式为f\left( x\right) =x^{k}，是一个k次多项式。那么原式就为k+1次多项式。}\)
\(\\\)
\(\large{那么我们就是去求这个k+1次的多项式f(n)的值，但是插值的复杂度是\text{O}(\text{k}^2)，显然会T，但是这个题目我们可以随便取k+2个点进行插值，那么我们当然去找有特殊性质的点来考虑优化。\\于是选择从1起连续的k+2个点，那么原式子可以化为\\ f\left( n\right) =\sum \limits^{k+2}_{i=1}\left( \sum \limits^{i}_{j=1}j^{k}\right) \times \dfrac {pre_{i-1} \times suf_{i+1}}{fac_i \times fac_{k + 2 - i}} \\其中pre_i表示k-i的前缀积，suf_i表示k到i的后缀积，fac_i表示i!。那么这样可以现预处理出pre、suf与fac，然后y_i每次logk维护即可，循环k次，故时间复杂度为\text{O}(\text{klogk})。}\)
\(\\\)
\(\Large\textbf{Code: }\)

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

const int N = 1e6 + 5;
const int p = 1e9 + 7;

int n, k, suf[N], pre[N], fac[N];

void read(int &x) {
	int flg = 1; x = 0;
	char c = getchar();
	while (!isdigit(c)) { if (c == '-') flg = -1; c = getchar(); }
	while (isdigit(c)) x = x * 10 + c - '0', c = getchar();
	x *= flg;
}

int pow(int x, int y) {
	ll a = x, ret = 1;
	for (; y ; y >>= 1, a = a * a % p) if (y & 1) ret = ret * a % p;
	return ret;
}

int main() {
	read(n), read(k);
	suf[k + 3] = pre[0] = fac[0] = 1;
	for (int i = 1; i <= k + 2; ++i) pre[i] = 1ll * pre[i - 1] * (n - i) % p;
	for (int i = k + 2; i >= 1; --i) suf[i] = 1ll * suf[i + 1] * (n - i) % p;
	for (int i = 1;  i<= k + 2; ++i) fac[i] = 1ll * fac[i - 1] * i % p;
	int y = 0, f, s, ans = 0;
	for (int i = 1; i <= k + 2; ++i) {
		y = (1ll * y + pow(i, k)) % p;
		s = 1ll * pre[i - 1] * suf[i + 1] % p;
		f = 1ll * fac[i - 1] * ((k - i) & 1 ? -1ll : 1ll) * fac[k + 2 - i] % p;
		ans = (1ll * ans + 1ll * y * s % p * pow(f, p - 2) % p + p) % p;
	}
	printf("%d\n", ans);
	return 0;
}

\(\large{至于手推\sum \limits^{n}_{i=1}i^{k}的公式，我还没学，咕咕咕。}\)