拉格朗日插值学习笔记

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

(\)
(Large extbf{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
", ans);
	return 0;
} 

(\)
(large{再说一下拓展，如何运用拉格朗提插值在 ext{O}( ext{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} 是一个 ext {p} 阶等差数列的充要条件是数列的通项 a_n 为 ext {n} 的一个 ext {p} 次多项式。\证明：\有了上面那个定理，我们只关心最高项即可。\设数列 { a_n} 的通项 a_n是一个关于 ext {n} 的 ext {v} 次多项式，即fleft( x ight) =sum limits^{ u }_{i=0}u_{i}cdot x^{i}。\令数列b {n}为a {n}的一阶差分，那么b {n}的通项公式为fleft( x+1 ight) -fleft( x ight) =sum limits^{v}_{i=0}u_{i} imes left( x+1 ight) ^{i}-sum limits^{v}_{i=0}u_{i} imes x_{i}。\因为只考虑最高项，所以式子为u_v imes (x+1)^v - u_v imes x^v，然后用二项式定理把(x+1)^v展开，只考虑最高项，显然系数为1，那么式子就变成了u_v imes x^v - u_v imes 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}, sumlimits^{2}_{i=1}i^{k}... sumlimits^{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，显然这个数列的通项公式为fleft( x ight) =x^{k}，是一个k次多项式。那么原式就为k+1次多项式。})
(\)
(large{那么我们就是去求这个k+1次的多项式f(n)的值，但是插值的复杂度是 ext{O}( ext{k}^2)，显然会T，但是这个题目我们可以随便取k+2个点进行插值，那么我们当然去找有特殊性质的点来考虑优化。\于是选择从1起连续的k+2个点，那么原式子可以化为\ fleft( n ight) =sum limits^{k+2}_{i=1}left( sum limits^{i}_{j=1}j^{k} ight) imes dfrac {pre_{i-1} imes suf_{i+1}}{fac_i imes fac_{k + 2 - i}} \其中pre_i表示k-i的前缀积，suf_i表示k到i的后缀积，fac_i表示i!。那么这样可以现预处理出pre、suf与fac，然后y_i每次logk维护即可，循环k次，故时间复杂度为 ext{O}( ext{klogk})。})
(\)
(Large extbf{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
", ans);
	return 0;
}

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