首页 > 其他 > 详细

拉格朗日插值学习笔记

时间:2020-04-15 01:27:39      阅读:90      评论:0      收藏:0      [点我收藏+]

\(\large{对于一个关于x的n次多项式f(x),若知道其中的n+1个点,拉格朗日插值可以在\text{O}(\text{n}^2)的时间复杂度内求出f(k)。}\)
\(\\\)
\(\large{式子是这样的:\\ f\left( k\right) =\sum \limits^{n}_{i=0}y_{i_{j}}\prod \limits^{n}_{j=0,j\neq i}\dfrac {k-x_{j}}{x_{i}-x_{j}}\\关于正确性的证明,我还没理解。贴一下大佬的说法。}\)
\(\\\)
技术分享图片
\(\\\)
\(\Large\textbf{Code: }\)

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

const int N = 2e3 + 5;
const int p = 998244353;

int n, k, x[N], y[N];

void read(int &x) {
	int flg = 1; x = 0;
	char c = getchar();
	while (!isdigit(c)) { if (c == ‘-‘) flg = -1; c = getchar(); }
	while (isdigit(c)) x = x * 10 + c - ‘0‘, c = getchar();
	x *= flg;
}

int pow(int x, int y) {
	ll a = x, ret = 1;
	for (; y ; y >>= 1, a = a * a % p) if (y & 1) ret = ret * a % p;
	return  ret;
}

int main() {
	read(n), read(k);
	for (int i = 1; i <= n; ++i) read(x[i]), read(y[i]);
	int s, f, ans = 0;
	for (int i = 1; i <= n; ++i) {
		f = y[i], s = 1;
		for (int j = 1; j <= n; ++j) {
			if (i == j) continue;
			f = 1ll * f * ((k - x[j] + p) % p) % p, s = 1ll * s * ((x[i] - x[j] + p) % p) % p;
		}
		s = pow(s, p - 2);
		ans = (ans + 1ll * s * f % p) % p;
	}
	printf("%d\n", ans);
	return 0;
} 

\(\\\)
\(\large{再说一下拓展,如何运用拉格朗提插值在\text{O}(\text{klogk})的时间内求出\sum ^{n}_{i=1}i^{k},其中 n \leq 1e9。}\)
\(\\\)
\(\large{首先证明一下\sum ^{n}_{i=1}i^{k}这个式子是一个k+1次的多项式。太多了,待补,咕咕咕。}\)
\(\\\)
\(\large{那么我们就是去求这个k+1次的多项式f(n)的值,但是插值的复杂度是\text{O}(\text{n}^2),显然会T,但是这个题目我们可以随便取k+2个点进行插值,那么我们当然去找有特殊性质的点来考虑优化。\\于是选择从1起连续的k+2个点,那么原式子可以化为\\ f\left( n\right) =\sum \limits^{k+2}_{i=1}\left( \sum ^{i}_{j=1}j^{k}\right) \times \dfrac {pre_{i-1} \times suf_{i+1}}{fac_i \times fac_{k + 2 - i}} \\其中pre_i表示k-i的前缀积,suf_i表示k到i的后缀积,fac_i表示i!。那么这样可以现预处理出pre、suf与fac,然后y_i每次logk维护即可,循环k次,故时间复杂度为\text{O}(\text{klogk})。}\)
\(\\\)
\(\Large\textbf{Code: }\)

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

const int N = 1e6 + 5;
const int p = 1e9 + 7;

int n, k, suf[N], pre[N], fac[N];

void read(int &x) {
	int flg = 1; x = 0;
	char c = getchar();
	while (!isdigit(c)) { if (c == ‘-‘) flg = -1; c = getchar(); }
	while (isdigit(c)) x = x * 10 + c - ‘0‘, c = getchar();
	x *= flg;
}

int pow(int x, int y) {
	ll a = x, ret = 1;
	for (; y ; y >>= 1, a = a * a % p) if (y & 1) ret = ret * a % p;
	return ret;
}

int main() {
	read(n), read(k);
	suf[k + 3] = pre[0] = fac[0] = 1;
	for (int i = 1; i <= k + 2; ++i) pre[i] = 1ll * pre[i - 1] * (n - i) % p;
	for (int i = k + 2; i >= 1; --i) suf[i] = 1ll * suf[i + 1] * (n - i) % p;
	for (int i = 1;  i<= k + 2; ++i) fac[i] = 1ll * fac[i - 1] * i % p;
	int y = 0, f, s, ans = 0;
	for (int i = 1; i <= k + 2; ++i) {
		y = (1ll * y + pow(i, k)) % p;
		s = 1ll * pre[i - 1] * suf[i + 1] % p;
		f = 1ll * fac[i - 1] * ((k - i) & 1 ? -1ll : 1ll) * fac[k + 2 - i] % p;
		ans = (1ll * ans + 1ll * y * s % p * pow(f, p - 2) % p + p) % p;
	}
	printf("%d\n", ans);
	return 0;
}

\(\large{至于手推\sum \limits^{n}_{i=1}i^{k}的公式,我还没学,咕咕咕。}\)

拉格朗日插值学习笔记

原文:https://www.cnblogs.com/Miraclys/p/12702236.html

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!