更好实现的简单插值法—

### 一、拉格朗日插值
*所有数列都有规律——拉格朗日*

比起牛顿插值，拉格朗日插值的名号要响亮的多。
而这个插值的想法和做法也十分简单：

设有$n$个点$(x_i,y_i)$，拉格朗日告诉我们一定有一个最多$n-1$次的多项式穿过这些点。
$$\sum_{p}y_p\prod_{i\neq p}\frac{x-x_i}{x_p-x_i}$$
多么简单的道理！对于每一个$p$来说，只要$x$取到其它点的$x_i$那么整个连乘积为$0$，就没有了这一项$y_p$，只有在$x$取到$x_p$的时候后面的系数为$1$。

我们可以把每一个二项式展开，把$n$个二项式加起来，就可以求得所求多项式。总复杂度$O(n^3)$
但是通常情况我们只是求单点的值
直接带到没有化简的式子中去，加上预处理就可以$O(n^2)$

现在我们看到了这个算法的瓶颈——虽然单点计算很迅速，但是要求得目标多项式的话代码复杂且时间较劣。

因此，这篇文章将介绍一个实现简单且可以计算出原多项式的插值算法。

### 二、牛顿插值

牛顿插值是一个增量算法。
牛顿的思路就比较不一样了，假设已知前$n-1$个点的多项式$f_{n-1}(x)$，是否可以求出前$n$个点的多项式$f_n(x)$呢？

首先引入**差商的概念**
定义一个点集的差商如下：
$0$阶差商$f[x_0]=f(x_0)$
$1$阶差商$f[x_0,x_1]=\frac{f[x_1]-f[x_0]}{x_1-x_0}$
$2$阶差商$f[x_0,x_1,x_2]=\frac{f[x_1,x_2]-f[x_0,x_1]}{x_2-x_0}$
$....$
$n$阶差商$f[x_0,x_1,...,x_n]=\frac{f[x_1,x_2,...,x_n]-f[x_0,x_1,...,x_{n-1}]}{x_n-x_0}$

这个东西到底有什么用呢

我们考虑存在未知点$x$的差商，简单移项可得：
$f[x]=f(x)$
$f[x,x_0]({x-x_0})+f[x_0]={f[x]}$
$f[x,x_0,x_1]({x-x_1})+f[x_0,x_1]={f[x,x_0]}$
$....$
将后面的项源源不断地带入前面：
$f(x)=f[x_0]+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+....+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})+f[x,x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_n)$

因为最后一项无论$x$取哪一个都会是$0$，可以去掉。

最后就是
$f(x)=f[x_0]+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+....+f[x_0,x_1,...,x_n](x-x_0)(x-x_1)...(x-x_{n-1})$

这就是牛顿插值公式。
实现的时候便可以增量构造，我们时刻维护$x_n$结尾的$n$个差商和多项式$\prod(x-x_i)$即可，差商用定义的公式递推，多项式因为每次乘二项式，直接递推，两个东西都可以$O(n)$维护。
加上求逆元(因为值域较大无法预处理)，复杂度为$O(n^2logn)$

代码实现如下：

```
#include<cstdio>
#define LL long long
using namespace std;
const int mod = 998244353;
const int maxn = 1e5 + 5;
int dt[2][maxn], px[maxn], py[maxn], pcnt, now;
int ans[maxn], plg[maxn], inv[maxn];

inline int mul(const int &a, const int &b) {
	return 1ll * a * b % mod;
}

inline int dec(const int &a, const int &b) {
	return a - b < 0 ? a - b + mod : a - b;
}

inline int add(const int &a, const int &b) {
	return a + b >= mod ? a + b - mod : a + b;
}

inline void Add(int &x, const int &a) {
	x += a; if(x >= mod) x -= mod;
}

LL fpow(LL a, int b) {
	LL ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		a = a * a % mod, b >>= 1;
	}
	return ans;
}

void init(int x, int y) {
	now = 0;
	dt[now][0] = y;
	px[++pcnt] = x;
	py[pcnt] = y, plg[0] = 1;
	ans[0] = y;
}

void insert(int x, int y) {
	px[++pcnt] = x;
	py[pcnt] = y, now ^= 1, dt[now][0] = y;
	//维护差商 
	for(register int i = 1; i < pcnt; ++i) {
		dt[now][i] = mul(dec(dt[now][i - 1], dt[now ^ 1][i - 1]), 
						fpow(dec(px[pcnt], px[pcnt - i]), mod - 2));
	}
	//维护多项式 
	for(register int i = pcnt - 1; i >= 1; --i) {
		plg[i] = add(plg[i - 1], mul((mod - px[pcnt - 1]), plg[i]));
	}
	plg[0] = mul(mod - px[pcnt - 1], plg[0]);
	//累计答案 
	for(register int i = 0; i < pcnt; ++i)
		Add(ans[i], mul(plg[i], dt[now][pcnt - 1]));
}

int main() {
	int x, y, n,ret = 0;
	scanf("%d", &n), inv[0] = 1;
	scanf("%d%d", &x, &y);
	for(register int i = 1; i <= n; ++i)
		inv[i] = (mod - mod / i) * inv[mod % i];
	init(x, y);
	for(register int i = 2; i <= n; ++i)
		scanf("%d%d", &x, &y), insert(x, y);
	for(register int i = 0; i < pcnt; ++i) {
		printf("%d ", ans[i]);
	}
	return 0;
} 
```

Rockdu's Blog

导航

最近发表

大佬们的友链