Lagrange 插值法学习笔记

插值

定义

给定 \(n\) 个离散数据点（称为节点）\((x_k,y_k),k=1,2,...,n\)。对于 \(x,(x \neq x_k,k=1,2,...,n)\) ，求 \(x\) 所对应的 \(y\) 的值称为内插。

\(f(x)\) 为定义在区间 \([a,b]\) 上的函数。\(x_1,x_2,x_3,\ldots,x_n\) 为 \([a,b]\) 上 \(n\) 个互不相同的点， \(G\) 为给定的某一函数类。若 \(G\) 上有函数 \(g(x)\) 满足：

\[g(x_i) = f(x_i),k =1,2,\ldots,n \]
则称 \(g(x)\) 为 \(f(x)\) 关于节点 \(x_1,x_2,x_3,\ldots,x_n\) 在 \(G\) 上的插值函数。

—— 插值（离散数学名词）_百度百科 (baidu.com)

应用

插值：用来填充图像变换时像素之间的空隙。

—— 插值（离散数学名词）_百度百科 (baidu.com)

基本类型

多项式
埃尔米特
分段
三角函数

Lagrange 插值法

用于求解多项式插值，在 OI 中具有广泛运用，经常用于拟合 \(n\) 次多项式。

分析

设 \(n\) 个离散数据点为 \(P_1(x_1,y_1),P_2(x_2,y_2),\ldots,P_n(x_n,y_n)\)，要求构造的函数为 \(f(x)\)。

构造 \(n\) 个函数 \(f_1(x),f_2(x),\dots,f_n(x)\)，对于其中第 \(i\) 个函数 \(f_i(x)\)，其图像过：（\(P_j'(x_j,0)\) 代表第 \(i\) 个点在 \(x\) 轴上的投影）

\[\begin{cases} P_j'(x_j,0) , j \neq i \\ P_i(x_i,x_j) \end{cases} \]

那么 \(f(x)\) 就等于 \(\sum_{i=1}^n f_i(x)\)。

可以简单地设：

\[f_i(x) = a_i \cdot \prod_{j=1}^n (x - x_j) [i \neq j] \]

将各个离散数据点代入即可得到：

\[a_i = \frac{y_i}{\prod_{j=1}^n (x_i - x_j) [i \neq j]} \\ \begin{aligned} f_i(x) & = a_i \cdot \prod_{j=1}^n (x - x_j) [i \neq j] \\ f_i(x) & = \frac{y_i}{\prod_{j=1}^n (x_i - x_j) [i \neq j]} \cdot \prod_{j=1}^n (x - x_j) [i \neq j] \\ f_i(x) & = y_i \prod_{j=1}^n \frac{x - x_j}{x_i - x_j} [i \neq j] \\ \end{aligned} \]

那么最后累加得到：

\[\begin{aligned} f(x) & = \sum_{i=1}^n f_i(x) \\ f(x) & = \sum_{i=1}^n y_i \prod_{j=1}^n \frac{x - x_j}{x_i - x_j} [i \neq j] \\ \end{aligned} \]

可以在 \(O(n^2)\) 的时间复杂度内实现 Lagrange 插值。

特殊情况：横坐标是连续整数的 Lagrange 插值

如果 \(n\) 个离散数据点的 \(x_1,x_2,x_3,\ldots,x_n\) 可以排成连续的整数，那么我们可以实现时间复杂度 \(O(n)\) 的插值。

设我们已知的离散数据点分别为 \(P_1(z,f(z)),P_2(z+1,f(z+1)),...,P_n(z+n-1,f(z+n-1))\)，那么代入 \(f(x)\) 得到：

\[\begin{aligned} f(x) & = \sum_{i=1}^n y_i \prod_{j=1}^n \frac{x - x_j}{x_i - x_j} [i \neq j] \\ f(x) & = \sum_{i=1}^n y_i \prod_{j=1}^n \frac{x - (z + j - 1)}{i-j} [i \neq j] \\ f(x) & = \sum_{i=1}^n y_i \frac{\frac{\prod_{j=1}^n[(x - z - 1) - j]}{(x - z - 1) - i}}{\prod_{j=1}^n(i-j)[i \neq j]} \\ f(x) & = \sum_{i=1}^n (-1)^{n-i} y_i \frac{{\prod_{j=1}^n[(x - z - 1) - j]}}{[(x - z - 1) - i](i-1)!(n-i)!} \\ \end{aligned} \]

那么只需要 \(O(n)\) 的预处理即可解决问题。

应用

【模板】拉格朗日插值 - 洛谷

这题就是直接给定离散数据点，并且要求插值。

#include<bits/stdc++.h>
#define INF 0x3f3f3f3f
#define Fi first
#define Se second
#define ll long long
#define Pii pair<int,int>
#define RCL(a,b,c,d) memset(a,b,sizeof(c)*(d))
#define FOR(i,a,b) for(int i(a);i<=(int)(b);++i)
#define DOR(i,a,b) for(int i(a);i>=(int)(b);--i)
#define tomax(a,...) ((a)=max({(a),__VA_ARGS__}))
#define tomin(a,...) ((a)=min({(a),__VA_ARGS__}))
#define EDGE(g,i,x,y) for(int i=(g).h[(x)],y=(g)[(i)].v;~i;y=(g)[i=(g)[i].nxt].v)
#define main Main();signed main(){ios::sync_with_stdio(0);cin.tie(0),cout.tie(0);return Main();}signed Main
using namespace std;
constexpr int N(2e3+10);
namespace Modular {
#define Mod 998244353
	template<class T1,class T2>constexpr auto add(T1 a,T2 b) {
		return a+b>=Mod?a+b-Mod:(a+b<0?a+b+Mod:a+b);
	}
	template<class T1,class T2>constexpr auto mul(T1 a,T2 b) {
		return (1ll*a*b%Mod+Mod)%Mod;
	}
	template<class T,class...Types>constexpr auto add(T a,Types...args) {
		return add(a,add(args...));
	}
	template<class T,class...Types>constexpr auto mul(T a,Types...args) {
		return mul(a,mul(args...));
	}
	template<class T1,class T2>T1 &toadd(T1 &a,T2 b) {
		return a=add(a,b);
	}
	template<class T1,class T2>T1 &tomul(T1 &a,T2 b) {
		return a=mul(a,b);
	}
	template<class T0,class T,class...Types>T0 &toadd(T0 &a,T b,Types...args) {
		return toadd(a,b),toadd(a,args...);
	}
	template<class T0,class T,class...Types>T0 &tomul(T0 &a,T b,Types...args) {
		return tomul(a,b),tomul(a,args...);
	}
	int Pow(int a,int b=Mod-2) {
		int res(1);
		for(a%=Mod;b;b>>=1,tomul(a,a))if(b&1)tomul(res,a);
		return res;
	}
} using namespace Modular;
int n,k;
struct Lagrange {
	vector<Pii> vec;
	void Init() {
		vec.clear();
	}
	void Scan(int n) {
		FOR(i,1,n) {
			int x,y;
			cin>>x>>y,vec.push_back(Pii(x,y));
		}
	}
	int Inter(int k) {
		int sum(0);
		for(const Pii &a:vec) {
			Pii cur(a.Se,1);
			for(const Pii &b:vec)if(a.Fi!=b.Fi)tomul(cur.Fi,add(k,Mod-b.Fi)),tomul(cur.Se,add(a.Fi,Mod-b.Fi));
			toadd(sum,mul(cur.Fi,Pow(cur.Se)));
		}
		return sum;
	}
} Lag;
signed main() {
	cin>>n>>k,Lag.Scan(n),cout<<Lag.Inter(k)<<endl;
	return 0;
}

The Sum of the k-th Powers - Codeforces

给出 \(n,k\)，求解 \(\sum_{i=1}^n i^k\)。

设数列 \(a_x = \sum_{i=1}^x i^k\)。

通过高阶等差的性质，我们可知当数列 \(b_i = a_{i+1} - a_{i}\) 的通项公式是一个 \(k\) 次多项式的时候，那么 \(\{a_i\}\) 的通项公式就是一个 \(k+1\) 次多项式，所以我们求出 \(a_1,a_2,a_3,\ldots,a_{k+2}\)，就可以进行 \(O(k)\) 的插值。

怎么求 \(a_1,a_2,a_3,\ldots,a_{k+2}\) 呢？

很简单，用欧拉线性筛求解 \(i^k,i=1,2,\ldots,k+2\)，再用前缀和累加，时间复杂度为 \(O(\frac{n\log_2{P}}{\ln{n}})\)，其中 \(P = 10^9 + 7\)，为题中模数。

#include<bits/stdc++.h>
#define INF 0x3f3f3f3f
#define ll long long
#define RCL(a,b,c,d) memset(a,b,sizeof(c)*(d))
#define FOR(i,a,b) for(int i(a);i<=(int)(b);++i)
#define DOR(i,a,b) for(int i(a);i>=(int)(b);--i)
#define tomax(a,...) ((a)=max({(a),__VA_ARGS__}))
#define tomin(a,...) ((a)=min({(a),__VA_ARGS__}))
#define EDGE(g,i,x,y) for(int i=(g).h[(x)],y=(g)[(i)].v;~i;y=(g)[i=(g)[i].nxt].v)
#define main Main();signed main(){ios::sync_with_stdio(0);cin.tie(0),cout.tie(0);return Main();}signed Main
using namespace std;
constexpr int N(1e6+10);
namespace Modular {
#define Mod 1000000007
	int inv[N],fac[N],ifac[N];
	template<class T1,class T2>constexpr auto add(T1 a,T2 b) {
		return a+b>=Mod?a+b-Mod:(a+b<0?a+b+Mod:a+b);
	}
	template<class T1,class T2>constexpr auto mul(T1 a,T2 b) {
		return (1ll*a*b%Mod+Mod)%Mod;
	}
	template<class T,class...Types>constexpr auto add(T a,Types...args) {
		return add(a,add(args...));
	}
	template<class T,class...Types>constexpr auto mul(T a,Types...args) {
		return mul(a,mul(args...));
	}
	template<class T1,class T2>T1 &toadd(T1 &a,T2 b) {
		return a=add(a,b);
	}
	template<class T1,class T2>T1 &tomul(T1 &a,T2 b) {
		return a=mul(a,b);
	}
	template<class T0,class T,class...Types>T0 &toadd(T0 &a,T b,Types...args) {
		return toadd(a,b),toadd(a,args...);
	}
	template<class T0,class T,class...Types>T0 &tomul(T0 &a,T b,Types...args) {
		return tomul(a,b),tomul(a,args...);
	}
	int Pow(int a,int b=Mod-2) {
		int res(1);
		for(a%=Mod; b; b>>=1,tomul(a,a))if(b&1)tomul(res,a);
		return res;
	}
	void Init(int lim=N-5) {
		fac[0]=1;
		FOR(i,1,lim)fac[i]=mul(fac[i-1],i);
		ifac[lim]=Pow(fac[lim]);
		DOR(i,lim,1)ifac[i-1]=mul(ifac[i],i);
		FOR(i,1,lim)inv[i]=mul(fac[i-1],ifac[i]);
	}
} using namespace Modular;
int n,k;
struct Lagrange {
	int y[N],pre[N],suf[N];
	vector<int> Pr;
	bitset<N> ip;
	void Eul(int lim,int k) {
		Pr.clear(),ip.reset(),ip[0]=ip[1]=1,y[1]=1;
		FOR(i,2,lim) {
			if(!ip[i])Pr.push_back(i),y[i]=Pow(i,k);
			for(const int &j:Pr) {
				if((ll)i*j>lim)break;
				ip[i*j]=1,y[i*j]=mul(y[i],y[j]);
				if(!(i%j))break;
			}
		}
		FOR(i,2,lim)toadd(y[i],y[i-1]);
	}
	int Inter(int x,int k) {
		Eul(k+2,k);
		if(x<=k+2)return y[x];
		int ans(0);
		pre[0]=suf[k+3]=1;
		FOR(i,1,k+2)pre[i]=mul(pre[i-1],x-i);
		DOR(i,k+2,1)suf[i]=mul(suf[i+1],x-i);
		FOR(i,1,k+2)toadd(ans,mul((k+2-i)&1?Mod-1:1,y[i],pre[i-1],suf[i+1],ifac[i-1],ifac[k+2-i]));
		return ans;
	}
} Lag;
signed main() {
	Init(),cin>>n>>k;
	cout<<Lag.Inter(n,k)<<endl;
	return 0;
}

[NOIP2020] 微信步数

//#define Plus_Cat "walk"
#include<bits/stdc++.h>
#define INF 0x3f3f3f3f
#define ll long long
#define RCL(a,b,c,d) memset(a,b,sizeof(c)*(d))
#define FOR(i,a,b) for(int i(a);i<=(int)(b);++i)
#define DOR(i,a,b) for(int i(a);i>=(int)(b);--i)
#define tomin(a,...) ((a)=min({(a),__VA_ARGS__}))
#define tomax(a,...) ((a)=max({(a),__VA_ARGS__}))
#define EDGE(g,i,x,y) for(int i(g.h[x]),y(g[i].v);~i;y=g[i=g[i].nxt].v)
using namespace std;
constexpr int N(5e5+10),M(20);
namespace Modular {
#define Mod 1000000007
	int fac[M],ifac[M];
	template<class T1,class T2>constexpr auto add(T1 a,T2 b) {
		return a+b>=Mod?a+b-Mod:(a+b<0?a+b+Mod:a+b);
	}
	template<class T1,class T2>constexpr auto mul(T1 a,T2 b) {
		return (1ll*a*b%Mod+Mod)%Mod;
	}
	template<class T,class...Types>constexpr auto add(T a,Types...args) {
		return add(a,add(args...));
	}
	template<class T,class...Types>constexpr auto mul(T a,Types...args) {
		return mul(a,mul(args...));
	}
	template<class T1,class T2>T1 &toadd(T1 &a,T2 b) {
		return a=add(a,b);
	}
	template<class T1,class T2>T1 &tomul(T1 &a,T2 b) {
		return a=mul(a,b);
	}
	template<class T0,class T,class...Types>T0 &toadd(T0 &a,T b,Types...args) {
		return toadd(a,b),toadd(a,args...);
	}
	template<class T0,class T,class...Types>T0 &tomul(T0 &a,T b,Types...args) {
		return tomul(a,b),tomul(a,args...);
	}
	int Pow(int a,int b=Mod-2) {
		int res(1);
		for(a%=Mod; b; b>>=1,tomul(a,a))if(b&1)tomul(res,a);
		return res;
	}
	void Init(int lim=M-5) {
		fac[0]=1;
		FOR(i,1,lim)fac[i]=mul(fac[i-1],i);
		ifac[lim]=Pow(fac[lim]);
		DOR(i,lim,1)ifac[i-1]=mul(ifac[i],i);
	}
} using namespace Modular;
int n,m,ans;
int a[M],b[M],c[N],d[N],e[M],h[M],w[M];
int f[2][M],l[N][M],r[N][M];
bool No_Solution() {
	static int l[M],r[M];
	RCL(l,0,l,1),RCL(r,0,r,1);
	FOR(i,1,n)e[c[i]]+=d[i],tomax(r[c[i]],e[c[i]]),tomin(l[c[i]],e[c[i]]);
	FOR(i,1,m)if(r[i]-l[i]>=w[i])return 0;
	FOR(i,1,m)if(e[i])return 0;
	return 1;
}
namespace Lagrange {
	int y[M],suf[M],pre[M];
	vector<int> Pr;
	bitset<M> ip;
	void Eul(int lim,int k) {
		Pr.clear(),ip.reset(),ip[0]=ip[1]=1,y[1]=1;
		FOR(i,2,lim) {
			if(!ip[i])Pr.push_back(i),y[i]=Pow(i,k);
			for(const int &j:Pr) {
				if(i*j>lim)break;
				ip[i*j]=1,y[i*j]=mul(y[i],y[j]);
				if(!(i%j))break;
			}
		}
		FOR(i,2,lim)toadd(y[i],y[i-1]);
	}
	//\sum_{i=1}^x i^k
	int F(int x,int k) {
		Eul(k+2,k);
		if(x<=k+2)return y[x];
		int ans(0);
		pre[0]=suf[k+3]=1;
		FOR(i,1,k+2)pre[i]=mul(pre[i-1],x-i);
		DOR(i,k+2,1)suf[i]=mul(suf[i+1],x-i);
		FOR(i,1,k+2)toadd(ans,mul((k+2-i)&1?Mod-1:1,y[i],pre[i-1],suf[i+1],ifac[i-1],ifac[k+2-i]));
		return ans;
	}
} using namespace Lagrange;
int main() {
#ifdef Plus_Cat
	freopen(Plus_Cat ".in","r",stdin),freopen(Plus_Cat ".out","w",stdout);
#endif
	scanf("%d%d",&n,&m),Init();
	FOR(i,1,m)scanf("%d",&w[i]);
	FOR(i,1,n)scanf("%d%d",&c[i],&d[i]);
	if(No_Solution())return puts("-1"),0;
	RCL(e,0,e,1);
	FOR(i,1,n) {
		e[c[i]]+=d[i];
		FOR(j,1,m)l[i][j]=l[i-1][j],r[i][j]=r[i-1][j];
		tomin(l[i][c[i]],e[c[i]]),tomax(r[i][c[i]],e[c[i]]);
	}
	FOR(j,1,m)a[j]=w[j]-(r[n][j]-l[n][j]);
	FOR(i,0,n) {
		int s(1);
		FOR(j,1,m)tomul(s,max(0,w[j]-(r[i][j]-l[i][j])));
		toadd(ans,s);
	}
	FOR(i,1,n)FOR(j,1,m)tomin(l[i][j]+=e[j]-l[n][j],0),tomax(r[i][j]+=e[j]-r[n][j],0);
	FOR(j,1,m)b[j]=r[n][j]-l[n][j];
	int last(-1);
	FOR(i,1,n) {
		RCL(f[0]+1,0,int,m),f[0][0]=1;
		int tmp(INF);
		FOR(j,1,m) {
			int K(-b[j]),B(a[j]-(r[i][j]-l[i][j]));
			if(B<=0)return printf("%d\n",ans),0;
			if(b[j]>0)tomin(tmp,B/b[j]);
			FOR(k,0,m) {
				f[j&1][k]=mul(f[!(j&1)][k],B);
				if(k)toadd(f[j&1][k],mul(f[!(j&1)][k-1],K));
			}
		}
		if(last!=tmp) {
			last=tmp,h[0]=tmp+1;
			FOR(j,1,m)h[j]=F(tmp,j);
		}
		FOR(j,0,m)toadd(ans,mul(h[j],f[m&1][j]));
	}
	printf("%d\n",ans);
	return 0;
}

posted @ 2024-10-16 15:52 Add_Catalyst 阅读(68) 评论(0) 收藏举报

刷新页面返回顶部

珈猫的小窝

欢迎来到珈猫的小窝

Lagrange 插值法学习笔记

Lagrange 插值法学习笔记

插值

定义

应用

基本类型

Lagrange 插值法

分析

特殊情况：横坐标是连续整数的 Lagrange 插值

应用

【模板】拉格朗日插值 - 洛谷

The Sum of the k-th Powers - Codeforces

[NOIP2020] 微信步数

公告

珈猫的小窝

欢迎来到珈猫的小窝

Lagrange 插值法 学习笔记

Lagrange 插值法 学习笔记

插值

定义

应用

基本类型

Lagrange 插值法

分析

特殊情况：横坐标是连续整数的 Lagrange 插值

应用

【模板】拉格朗日插值 - 洛谷

The Sum of the k-th Powers - Codeforces

[NOIP2020] 微信步数

公告

Lagrange 插值法学习笔记

Lagrange 插值法学习笔记