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

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

介绍

拉格朗日插值是用来求高次多项式的一种方法。我们知道假设已经给定了 \(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