拉格朗日插值法 以2019南昌icpc B-Polynomial

拉格朗日插值法

拉格朗日插值法

对任意一个$x$求$f(x)$

给定一个$n$次多项式（共有$n+1$项）
$$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$$

对于函数，若已知其任意$n+1$个x的取值以及其所对应的$f(x)$(即为$y$)，可以利用拉格朗日插值来求其他任意$X$的$f(X)$

$$f(X)=L(X)=\sum _{i=0}^n{y}_i\ast li(x)=\sum _{i=0}^n{y}_i\ast \prod _{j̸=j}^{j<=n}\frac{X-x_j}{x_i-x_j}$$

以$f(1)=1,f(2)=4,f(3)=9$为例
$$l_1(x)=\frac{x-2}{1-2}\ast \frac{x-3}{1-3}$$ $$l_2(x)=\frac{x-1}{2-1}\ast \frac{x-3}{2-3}$$$$l_3(x)=\frac{x-1}{3-1}\ast \frac{x-2}{3-2}$$
$$f(X)=f(1)\ast l_1(X)+f(2)\ast l_2(x)+f(3)\ast l_3(x)$$$$=1\ast \frac{x-2}{1-2}\ast \frac{x-3}{1-3}+4\ast \frac{x-1}{2-1}\ast \frac{x-3}{2-3}+9\ast \frac{x-1}{3-1}\ast \frac{x-2}{3-2}$$$$=x^2$$
带入公式对于任意一个x是$O(n)$的，多次询问即为$O(n^2)$

已知$f(x)$为连续的求$f(n)$

$$l_0(x)=\frac{(x-1)(x-2)\cdots (x-n)}{(0-1)(0-2)\cdots (0-n)}$$$$l_1(x)=\frac{(x-0)(x-2)\cdots (x-n)}{(1-0)(1-2)\cdots (0-n)}$$$$\cdots \cdots$$$$l_n(x)=\frac{(x-0)(x-1)\cdots (x-(n-1))}{(n-0)(n-1)\cdots (n-(n-1))}$$

分子：
令$$M=(x-0)(x-1)(x-2)\cdots (x-n)$$
$l_i(x)$的分子为$M/(x-i)$

分母
$$l_0(x)->(-1{)}^n\ast 0!\ast n!$$$$l_1(x)->(-1{)}^{n-1}\ast 1!\ast (n-1)!$$$$l_2(x)->(-1{)}^{n-2}\ast 2!\ast (n-2)!$$
分母为：$(-1{)}^{n-i}\ast i!\ast (n-i)!$

n！的逆元
先用费马小定理求出$n!$的逆元
$$\frac{1}{(n-1)!}=\frac{1}{n!}\ast n\mathrm{m}\mathrm{o}\mathrm{d}mod$$

Polynomial

题意

已知$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n^n$
你知道$f(0),f(1),f(2),\cdots ,f(n)$
求${\sum }_{i=L}^{R}f(i)$

分析

已知
$$f(0)=a_0+a_{1}0+a_2{0}^2+\cdots +a_n{0}^n$$$$f(1)=a_0+a_{1}1+a_2{1}^2+\cdots +a_n{1}^n$$$$f(2)=a_0+a_{1}2+a_2{2}^2+\cdots +a_n{2}^n$$$$\cdots \cdots$$$$f(n)=a_0+a_{1}n+a_2{n}^2+\cdots +a_n{n}^n$$

显然，对$1^n+2^n+3^n+\cdots +x^n=b_0{x}^{n+1}$

$$S(n)=\sum _{i=1}^{n}f(i)=b_0+b_{1}n+b_2{n}^2+\cdots +b_n{n}^n+b_{n+1}{n}^{n+1}$$

先用拉格朗日插值推出$a_{n+1}$，再推出$S(0),S(1),\cdots ,S(n),S(n+1)$
再用拉格朗日插值推出$S(R)-S(L-1)$

代码

预处理

fac[0] = 1;		//阶乘
for (int i = 1; i <= maxn; i++)
	fac[i] = fac[i] * i % mod;
inv[1] = 1;		//逆元
for (int i = 2; i <= maxn; i++) 
	inv[i] = (mod - mod / i) * (LL)1 * inv[mod % i] % mod;
fac_inv[maxn] = quickpow(fac_inv[maxn], mod - 2);	//n！的逆元
for (int i = maxn - 1; i >= 1; i--)
	fac_inv[i] = fac_inv[i + 1] * (i + 1) % mod;

拉格朗日插值$f(n+1)$

void Get_f_next() {		//计算a[n+1]
	f[n + 1] = 0;
	LL M = 1;		//处理分子
	for (int i = 0; i <= n; i++)
		M = M * (n + 1 - i) % mod;
	LL ans = 0;		//拉格朗日插值
	for (int i = 0; i <= n; i++) {
		ans = M * inv[n + 1 - i] % mod		//分子
			* fac_inv[n-i] % mod * fac_inv[i] % mod; //分母
		if ((n - i) & 1)f[n + 1] = (f[n + 1] - ans * f[i] % mod + mod) % mod;	//yi*li
		else f[n + 1] = (f[n + 1] + ans * f[i] % mod + mod) % mod;
	}
}

拉格朗日插值$S(r)$

LL Get_S(LL r) {
	if (r <= n + 1)
		return sum[r];
	LL M = 1;
	for (int i = 0; i <= n + 1; i++)	//处理分子
		M = M * LL(r - i) % mod;			//分母
	LL ans = 0, res = 0;
	for (int i = 0; i <= n + 1; i++) {
		ans = M * inv[r - i] % mod //分子
			* fac_inv[n + 1 - i] % mod * fac_inv[i] % mod;	//分母
		if ((n + 1 - i) & 1) res = (res - ans * sum[i] % mod + mod) % mod;	//yi*li
		else res = (res + ans * sum[i] % mod + mod) % mod;
	}
	return res;
}

完整代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
using namespace std;
#pragma warning (disable:4996)
typedef long long LL;
const int maxn = 1005;
const LL mod = 9999991;
LL fac[maxn + 5], inv[mod + 100], fac_inv[maxn + 5];
LL quickpow(LL m, LL p) {
	LL res = 1;
	while (p) {
		if (p & 1) 
			res = res * m % mod;
		m = m * m % mod;
		p >>= 1;
	}
	return res;
}
void init() {
	fac[0] = 1;
	for (int i = 1; i <= maxn; i++)
		fac[i] = fac[i - 1] * LL(i) % mod;
	fac_inv[maxn] = quickpow(fac[maxn], mod - 2);
	for (int i = maxn - 1; i >= 0; i--)
		fac_inv[i] = fac_inv[i + 1] * (LL(i) + 1) % mod;
	inv[0] = inv[1] = 1;
	for (int i = 2; i < mod+5; i++)
		inv[i] = (mod - mod / i) * (LL)1 * inv[mod % i] % mod;
}
LL t, n, m;
LL f[maxn], sum[maxn];
void Get_f_next() {		//计算a[n+1]
	f[n + 1] = 0;
	LL M = 1;		//处理分子
	for (int i = 0; i <= n; i++)
		M = M * (n + 1 - i) % mod;
	LL ans = 0;		//拉格朗日插值
	for (int i = 0; i <= n; i++) {
		ans = M * inv[n + 1 - i] % mod		//分子
			* fac_inv[n-i] % mod * fac_inv[i] % mod; //分母
		if ((n - i) & 1)f[n + 1] = (f[n + 1] - ans * f[i] % mod + mod) % mod;	//yi*li
		else f[n + 1] = (f[n + 1] + ans * f[i] % mod + mod) % mod;
	}
}
LL Get_S(LL r) {
	if (r <= n + 1)
		return sum[r];
	LL M = 1;
	for (int i = 0; i <= n + 1; i++)	//处理分子
		M = M * LL(r - i) % mod;			//分母
	LL ans = 0, res = 0;
	for (int i = 0; i <= n + 1; i++) {
		ans = M * inv[r - i] % mod //分子
			* fac_inv[n + 1 - i] % mod * fac_inv[i] % mod;	//分母
		if ((n + 1 - i) & 1) res = (res - ans * sum[i] % mod + mod) % mod;	//yi*li
		else res = (res + ans * sum[i] % mod + mod) % mod;
	}
	return res;
}
int main() {
	init();
	scanf("%lld", &t);
	while (t--) {
		scanf("%lld%lld", &n, &m);
		for (int i = 0; i <= n; i++) {
			cin >> f[i]; f[i] %= mod;
		}
		Get_f_next();
		sum[0] = f[0];
		for (int i = 1; i <= n + 1; i++)
			sum[i] = (sum[i - 1] + f[i]) % mod;
		LL l, r;
		while (m--) {
			scanf("%lld%lld", &l, &r);
			printf("%lld\n", (Get_S(r) - Get_S(l - 1) + mod) % mod);
		}
	}
}