拉格朗日插值学习笔记

最新推荐文章于 2024-08-26 23:12:24 发布

FSYo

最新推荐文章于 2024-08-26 23:12:24 发布

阅读量243

点赞数 1

分类专栏：拉格朗日插值

本文链接：https://blog.csdn.net/sslz_fsy/article/details/90549077

版权

拉格朗日插值专栏收录该内容

5 篇文章 0 订阅

订阅专栏

当给出了n个点的坐标，求经过这 n 个点的多项式，以及某个点的点值，可以在 O(N^2) 的时间内求出，而朴素的高斯消元需要O(N^3)

给出 (x_0, y_0),(x_1,y_1)...(x_n,y_n)

我们对于每一个 i，构造一个多项式 $l_i(x)=y_i\prod _{j\neq i}\frac{x-x_j}{x_i-x_j}$

发现，把 xi 带进去时，函数值是 yi, 而将其他的x带进去时，函数值都是0

于是我们把每个 i 的 li(x)相加，就是最后的多项式了

$f(x)=\sum_{i=0}^ny_i\prod _{j\neq i}\frac{x-x_j}{x_i-x_j}$

当要求 N 的函数值时，带进去就可以 n^2 做了

ll lagrange(ll k){
    ll ans = 0;
    for(int i=1; i<=n; i++){
        ll s1 = 1, s2 = 1;
        for(int j=1; j<=n; j++){
            if(i != j){
                s1 = mul(s1, k - x[j]);
                s2 = mul(s2, x[i] - x[j]);
            }
        }
        ans = add(ans, mul(y[i], mul(s1, power(s2, Mod-2))));
    } return ans;
}

如果x是连续的话，可以优化到O(n)

也就是求出所有(N - xj) 的积，然后不要的那项求一个逆元

分母是 (k-i)! *(i-1)!，要注意一下符号问题

ll lagrange(ll n, ll k){
	ll pre = 1, ans = 0;
	if(n <= k) return f[n];
	for(int i=1; i<=k; i++) pre = mul(pre, n - i);
	for(int i=1; i<=k; i++){
		ll flag = ((k-i) & 1) ? -1 : 1;
		ll inv1 = power(n - i, Mod-2);
		ll inv2 = (flag * power(mul(fac[i-1], fac[k-i]), Mod-2) + Mod) % Mod;
		ans = add(ans, mul(mul(inv1, inv2), mul(f[i], pre)));
	} return ans;
	
}

经典应用

首先讲一个经典应用：计算 $\sum_{i=1}^ni^k(n<=1e15)$

老祖宗告诉我们，这个东西是个以n为自变量的k + 1次多项式

然后直接带入k+1个点后用拉格朗日插值算即可，复杂度O(k)

也就是说，只要能判断要求的是一个多项式，我们就可以用拉格朗日来插

最后放两道道模板题 P4593 [TJOI2018]教科书般的亵渎

这道题主要就是求 $\sum_{i=1}^ni^k$ , 直接插就可以了

#include<bits/stdc++.h>
#define N 100
using namespace std;
const int Mod = 1e9+7;
typedef long long ll;
ll a[N], f[N], fac[N]; int T;
ll add(ll a, ll b){ return (a+b) % Mod;}
ll mul(ll a, ll b){ return (a*b) % Mod;}
ll dec(ll a, ll b){ return (a+Mod-b) % Mod;}
ll power(ll a, ll b){
	ll ans = 1; for(;b;b>>=1){
		if(b&1) ans = mul(ans, a);
		a = mul(a, a);
	} return ans;
}
ll lagrange(ll n, ll k){
	ll pre = 1, ans = 0;
	if(n <= k) return f[n];
	for(int i=1; i<=k; i++) pre = mul(pre, n - i);
	for(int i=1; i<=k; i++){
		ll flag = ((k-i) & 1) ? -1 : 1;
		ll inv1 = power(n - i, Mod-2);
		ll inv2 = (flag * power(mul(fac[i-1], fac[k-i]), Mod-2) + Mod) % Mod;
		ans = add(ans, mul(mul(inv1, inv2), mul(f[i], pre)));
	} return ans;
	
}
int main(){
	scanf("%d", &T);
	fac[1] = fac[0] = 1;
	for(int i=2; i<=N-5; i++) fac[i] = mul(fac[i-1], i);
	while(T--){
		ll n, m; scanf("%lld%lld", &n, &m);
		for(int i=1; i<=m; i++) scanf("%lld", &a[i]);
		sort(a+1, a+m+1);
		for(int i=1; i<=m+3; i++) f[i] = add(f[i-1], power(i, m+1));
		ll ans = 0;
		for(int i=0; i<=m; i++){
			ans = add(ans, lagrange(n-a[i], m+3));
			for(int j=i+1; j<=m; j++) ans = dec(ans, power(a[j]-a[i], m+1));
		} printf("%lld\n", (ans % Mod + Mod) % Mod);
	} return 0;
}

[tyvj1858] xlkxc

求 $\sum_{i=1}^n\sum_{j=1}^{a+id}\sum_{l=1}^j l^k$

首先求 $f(x)=\sum_{i=1}^xi^k,g(x)=\sum_{i=1}^xf(i),F(n)=\sum_{i=1}^ng(a+id)$

可以知道 f 是k+1项的多项式, g差分一下就是 f, 所以 g 是 k+2 项的多项式, F是k+3项的多项式

于是把 F的k+4项用 g 插出来, 然后暴力插 F的第n项

#include<bits/stdc++.h>
#define N 200
using namespace std;
typedef long long ll;
const ll Mod = 1234567891;
int T, k, a, n, d;
ll f[N], g[N], fac[N];
ll add(ll a, ll b){ return (a+b) % Mod;}
ll mul(ll a, ll b){ return (a*b) % Mod;}
ll power(ll a, ll b){
	ll ans = 1; for(;b;b>>=1){
		if(b&1) ans = mul(ans, a);
		a = mul(a, a);
	} return ans;
}
ll lagrange(ll *a, ll k, ll n){
	ll pre = 1, ans = 0;
	if(n <= k) return a[n];
	for(int i=1; i<=k; i++) pre = mul(pre, n - i);
	for(int i=1; i<=k; i++){
		ll inv1 = power(n - i, Mod-2);
		ll inv2 = power(mul(fac[i-1], fac[k-i]), Mod-2);
		ll flag = ((k-i) & 1) ? -1 : 1;
		ans = add(ans, flag * mul(mul(inv1, inv2), mul(a[i], pre)));
		ans = (ans % Mod + Mod) % Mod;
	} return ans;
}
int main(){
	scanf("%d", &T); fac[1] = fac[0] = 1; 
	for(int i=2; i<=N-5; i++) fac[i] = mul(fac[i-1], i);
	while(T--){
		scanf("%d%d%d%d", &k, &a, &n, &d);
		for(int i=1; i<=k+4; i++) f[i] = add(f[i-1], power(i, k));
		for(int i=1; i<=k+4; i++) f[i] = add(f[i], f[i-1]);
		for(int i=0; i<=k+4; i++) g[i] = lagrange(f, k+4, add(a, mul(i, d)));
		for(int i=1; i<=k+4; i++) g[i] = add(g[i], g[i-1]);
		printf("%lld\n", lagrange(g, k+4, n));
	} return 0;
}