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

前记

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

公式

\[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;
}