[学习笔记]拉格朗日插值

前记

这是一个很早就推出柿子但是一直没有做的一个知识。

公式

$$f(k)=\sum\limits_{i=1}^ny_i\prod\limits_{i\not=j}\frac{k-x_j}{x_i-x_j}$$

我有两种办法推导这个柿子，但是可惜这里空间太小写不下。

用法

计算 $\sum\limits_{i=1}^{n}i^k$

我们尝试从小的情况入手。

$\sum\limits_{i=1}^{n}i\ \ =\ \ \displaystyle\frac{n(n+1)}{2}$
$\sum\limits_{i=1}^{n}i ^ 2 \ \ = \ \ \displaystyle\frac{n(n+1)(2n+1)}{6}$

我们猜想，$\sum\limits_{i=1}^{n}i^k$ 是一个 $k+1$ 次多项式。很好的是，这是对的。

那么我们要计算 $\sum\limits_{i=1}^{n}i^k$ 相当于计算当 $x=n$ 时，一个多项式的值。这就是显然的拉格朗日插值了。可是如果用普通的拉格朗日插值，时间复杂度是 $\mathcal{O}(k^2)$ 的，这显然不可以接受，解决的办法也很简单，点值是随便我们带的，我们选择带入 $1,2,...,k+2$试试会不会有什么性质。

$f(n)=\sum\limits_{i=1}^{k+2}y_i\prod\limits_{i\not=j}\displaystyle\frac{k-j}{i-j}$

好像有什么东西在里面了。我们分子分母分开看。

对于分子，我们就是要预处理出 $\prod_{i\not=j}k-j$ 很好办，这类似于什么。

因为 $k$ 不会变把，所以这是一个固定序列挖掉一个格子求积。我们显然可以通过预处理前缀积和后缀积解决这个问题。

对于分母，依然是非常的好搞。

我们就是要处理 $\prod\limits_{i\not=j}\displaystyle\frac{1}{i-j}$，我们思考一下这个是什么。

好像会有负数，这对不对呢，是对的，因为原来我们也不需要保证 $x_i-x_j$ 是正数。

好像是一段连续的区间刚好挖掉了 $0$ 然后求积。

这段区间的范围就是 $(0,i-1]$ 另一段是 $[i-k, 0)$ 我们直接预处理阶乘。后面那个把符号全部拿出来，奇偶讨论一下就行了。

代码

#include <bits/stdc++.h>
#define int long long
using namespace std;
template <typename T>inline void read(T& t){t=0; register char ch=getchar(); register int fflag=1;while(!('0'<=ch&&ch<='9')) {if(ch=='-') fflag=-1;ch=getchar();}while(('0'<=ch&&ch<='9')){t=t*10+ch-'0'; ch=getchar();} t*=fflag;}
template <typename T,typename... Args> inline void read(T& t, Args&... args) {read(t);read(args...);}
const int N = 2e6 + 10, inf = 0x3f3f3f3f;

typedef long long ll;
int n, m;
const int P = 1e9 + 7;
ll pre[N], suf[N], x[N], y[N], fac[N], ifac[N], ans;

ll qpow(ll x, int y) {
	ll res = 1;
	while(y) {
		if(y & 1) res = res * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return res;
}

signed main() {
	read(n, m);
	ll res = 0;
	for(int i = 1; i <= m + 2; ++i) {
		res += qpow(i, m);
		res %= P;
		y[i] = res;
	}
	fac[0] = 1; for(int i = 1; i <= m + 10; ++i) fac[i] = fac[i - 1] * i % P;
	ifac[m + 10] = qpow(fac[m + 10], P - 2);
	for(int i = m + 9; i >= 0; --i) ifac[i] = ifac[i + 1] * (i + 1) % P;
	pre[0] = 1; for(int i = 1; i <= m + 2; ++i) pre[i] = pre[i - 1] * (n - i) % P;
	suf[m + 3] = 1; for(int i = m + 2; i >= 0; --i) suf[i] = suf[i + 1] * (n - i) % P;
	for(int i = 1; i <= m + 2; ++i) {
		// cout << '*' << m << endl;
		// cout << pre[i - 1] << ' ' << suf[i + 1] << ' ' << y[i] << ' ' << ifac[i - 1] << ' ' << ifac[m - i] << endl;
		// cout << pre[i - 1] * suf[i + 1] % P * y[i] % P * ifac[i - 1] % P * ifac[m - i] % P << endl;
		if((m - i) & 1) ans = (ans - pre[i - 1] * suf[i + 1] % P * y[i] % P * ifac[i - 1] % P * ifac[m + 2 - i] % P + P) % P;
		else ans = (ans + pre[i - 1] * suf[i + 1] % P * y[i] % P * ifac[i - 1] % P * ifac[m + 2 - i] % P) % P;
	}
	cout << ans << endl;
	return 0;
}