【学习笔记】拉格朗日插值

学习多项式的第一步。

参考资料：

attack的luogu博客

oi wiki拉格朗日插值

Apocryphal的luogu博客

1.拉格朗日插值法的简介

问题：

luogu P4781 【模板】拉格朗日插值

解法1：高斯消元

显然\(\deg \geqslant n\)的多项式有无穷个，因为根据高斯消元，一定会出现自由元。

直接把这\(n\)个点列出方程组，用高斯消元求解。

求出多项式后再求出\(f(k)\)即可。

时间复杂度\(O(n^3).\)

解法2：拉格朗日插值法

给出一个关于点\((x_i,y_i)\)的基函数：

\[g(k)=\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-x_j}{x_i-x_j} \]

容易发现：\(\forall j\not=i,g(x_j)=0.\)因为累乘中总有\(k=x_j\),使得\(k-x_j=x_j-x_j=0,g(x_j)=0.\)

对于\(j=i,g(x_j)=1.\)因为累乘中每一项均为\(\dfrac{x_i-x_j}{x_i-x_j}=1,g(x_j)=1.\)

于是我们的多项式就可以表示为：

\[f(k)=\sum_{i=1}^{n}y_i\times\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-x_j}{x_i-x_j} \]

这样\(\forall(x_i,y_i),f(x_i)=y_i\),也可以据此求出\(f(k).\)

大概因为要求逆元，所以时间复杂度为\(O(n^2+n\log n)=O(n^2).\)

\(\tt{code:}\)

#include <bits/stdc++.h>
#define ll long long
using namespace std;
int n;
const ll mod = 998244353;
ll x[2011], y[2011], k, ans;
ll ksm(ll s1, ll s2) {
	if(!s2) return 1;
	if(s2 & 1) return s1 * ksm(s1, s2 - 1) % mod;
	ll ret = ksm(s1, s2 / 2);
	return ret * ret % mod;
}
int main() {
	scanf("%d%lld", &n, &k);
	for(int i = 1; i <= n; i++) scanf("%lld%lld", &x[i], &y[i]);
	for(int i = 1; i <= n; i++) {
		ll sum = 1, mul = 1;
		for(int j = 1; j <= n; j++) {
			if(i == j) continue;
			sum *= (k - x[j]); sum %= mod; sum += mod; sum %= mod;
			mul *= (x[i] - x[j]); mul %= mod; mul += mod; mul %= mod;
		}
		sum *= ksm(mul, mod-2); sum %= mod;
		ans += (sum * y[i]) % mod; ans %= mod;
	}
	printf("%lld\n", ans);
	return 0; 
}

2.拉格朗日插值法的拓展

问题：

同上，加入\(x_i\)的值连续的限制。

解法：

假设\(x_i\in[1,n],\)则公式化为：

\[f(k)=\sum_{i=1}^{n}y_i\times\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-j}{i-j} \]

可以维护\(num=\prod_{1\leqslant j\leqslant n}k-j,\)于是分子部分化为\(\dfrac{num}{k-i},\)因为这么做要求逆元，带个log，很不优，所以可以预处理前缀和\(pre_i\)和后缀和\(suf_i.\)

对于分母，可以化为\(1\times 2\times···\times i-1\times (-1)\times (-2)\times ···\times(i-n)\)的形式。

实质上是一个阶乘，|分母|\(=fac_{i-1}\times fac_{n-i},\)然后再处理一下正负即可。

这样再加上一些\(O(n)\)的预处理\(pre,suf,fac,inv...\)，我们只需要\(O(n)\)即可完成。

3.拉格朗日插值法的具体应用&例题

【例1】CF622F The Sum of the k-th Powers

拉格朗日插值的经典模型，即求\(\sum_{i=1}^n i^k\).

这是个\(k+1\)次的多项式，感性证明一下：

考虑\(n\)次多项式\(g(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_0\).

它的一阶差分\(\Delta g(x)=g(x+1)-g(x)=a_n(x+1)^n-a_nx^n+...\).

通过二项式定理，我们发现这时的\(n\)次项会被消掉，即：

每做一次差分，多项式的度便会减一。

令\(f(x)=\sum_{i=1}^ni^k,\)则\(\Delta f(x)=f(x+1)-f(x)=(x+1)^k,\)为\(k\)次多项式。

所以\(f(x)\)为\(k+1\)次多项式。

至此，我们只需确定\(k+2\)个点，做一个插值，就可以求出\(f(n)\)了。

因为\(i\)连续，所以可以用到2.2中的拉格朗日插值拓展做到\(O(k\log k).\)

\(i^k\)也可以线性求出，这里不多作说明。

\(\tt{code:}\)

#include <bits/stdc++.h>

#define ll long long

using namespace std;

const int N = 1e6;
const ll mod = 1e9 + 7;

ll n, k, fac[N+11], pre[N+11], suf[N+11], ans;

ll ksm(ll s1, ll s2) {
	if(!s2) return 1ll;
	if(s2 % 2) return s1 * ksm(s1, s2-1) % mod;
	ll ret = ksm(s1, s2/2);
	return ret * ret % mod;
}

int main() {
	scanf("%lld%lld", &n, &k);
	ll sum = 0; 
	fac[0] = 1;
	for(int i = 1; i <= k+2; i++) 
		fac[i] = fac[i-1] * (ll)i % mod;
	pre[0] = suf[k+3] = 1;
	for(int i = 1; i <= k+2; i++)
		pre[i] = pre[i-1] * (n - i) % mod;
	for(int i = k+2; i >= 1; i--)
		suf[i] = suf[i+1] * (n - i) % mod;
	for(int i = 1; i <= k+2; i++) {
		sum += ksm(i, k); sum %= mod;
		ll num = pre[i-1] * suf[i+1] % mod;
		ll fm = fac[i-1] * fac[k+2-i] % mod;
		num *= sum; num %= mod;
		num *= ksm(fm, mod-2); num %= mod;
		if((k+2-i) % 2 == 0) ans += num;
		else ans -= num; 
		ans %= mod; ans = (ans + mod) % mod;
	}
	printf("%lld\n", ans);
	return 0;
}

【例2】luogu P4593 [TJOI2018]教科书般的亵渎

问题转化：先求出\(k,\)即“亵渎”的次数。忽略大于\(n\)的\(a_i\)（虽然不知道有没有），容易发现一次“亵渎”可以使一段血量连续的怪物全部死掉，所以\(k=m+1.\)

又因为怪物的血量互不相同，所以我们可以把求的值转化为\(\sum_{i=1}^n i^k,\)于是就和【例1】做法相似了。

注意如果有一段连续的\([x,n]\)的空位，不要计入答案。

代码就不放了。