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

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

介绍

拉格朗日插值是用来求高次多项式的一种方法。我们知道假设已经给定了 $n$ 个点，那么一定会有一个唯一确定的 $n-1$ 次的多项式，拉格朗日插值就是用来求这样的多项式的。

实现

给出拉格朗日插值的式子，假设已经给定了 $n$ 个点，分别为 $x_i,y_i,(1\le i \le n)$ 那么

$$f_x=\sum_{i=1}^ny_i\prod_{j=1,j\neq i}^n\frac{x-x_j}{x_i-x_j}$$

可以发现，对于第 $i$ 项，如果 $x$ 为 $x_i$ 则值为 $y_i$，否则为 $0$。因此把 $x_i$ 带入可得 $f_{}x_i=y_i$。暴力求解时间复杂度 $O(n^2)$

The Sum of the k-th Powers

思路

首先证明一个结论：对于一个通项公式为 $k$ 次的数列 $\{a_n\}$，每进行一次差分，生成新序列的通项公式的最高次数就会减 $1$。
首先设序列的通项公式为 $\sum_{i=1}^ku_i\times x^i$，其中 $u$ 为系数，则他的一阶差分的通项公式为

$$f_{x+1}-f_x=\sum_{i=1}^ku_i\times{((x+1)^i-x^i)}$$

将其用二项式定理拆开，然后发现 $x^k-x^k=0$ 所以最高次项约没了，因此得证。
回到题目，当前数列的通项公式为 $\sum_{i=1}^xi^k$，他的一阶差分的通项公式为

$$\sum_{i=1}^{x+1}i^k-\sum_{i=1}^xi^k=(x+1)^k$$

发现是一个 $k$ 次式，所以原始数列的通项公式为 $k+1$ 次。
所以我们可以先构造 $k+2$ 个点（直接暴力构造就行），然后使用拉格朗日插值。
考虑到基本公式 $f_x=\sum_{i=1}^ny_i\prod_{j=1,j\neq i}^n\frac{x-x_j}{x_i-x_j}$， 发现 $x_j,x_i$ 都是连续的，所以

• 分母上为 $x_i-1\times x_i-2... \times1\times-1...\times{x_i-(k+2)}$，发现为 $(x_i-1)!\times{(k+2-x_i)!}$
  注意需要判断符号如果 $(k-x_i+2)$ 如果为奇数需要再乘 $-1$。
• 分子用一个前后缀维护就行。

code

#include<bits/stdc++.h>
#define N 1000010
#define int long long
#define mod 1000000007
using namespace std;
int n, k, pre[N], suc[N], frac[N], y = 0, ans = 0;
int fast(int a, int p){
	int s = 1;
	while(p){
		if(p & 1) s = s * a % mod;
		a = a * a % mod;
		p = p >> 1;
	}
	//cout << s << endl;
	return s; 
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	cin >> n >> k;
	pre[0] = 1ll, suc[k + 3] = 1ll, frac[0] = 1ll;
	for(int i = 1; i <= k + 2; i++) frac[i] = frac[i - 1] * i % mod;
	for(int i = 1; i <= k + 2; i++) pre[i] = pre[i - 1] * (n - i) % mod;
	for(int i = k + 2; i >= 1; i--) suc[i] = suc[i + 1] * (n - i) % mod;
	for(int i = 1; i <= k + 2; i++){
		y = (y + fast(i, k)) % mod;
		int a = pre[i - 1] * suc[i + 1] % mod;
		int b = frac[i - 1] * frac[k + 2 - i] % mod  * ((k - i) & 1 ? -1ll : 1ll) % mod;
		ans = (ans + a % mod * fast(b, mod - 2) % mod * y % mod) % mod;
	}
	cout << (ans + mod) % mod;
	return 0;
}
``cpp