多项式插值法

最新推荐文章于 2024-08-01 11:03:33 发布

JK Chen

最新推荐文章于 2024-08-01 11:03:33 发布

阅读量2.2k

点赞数 1

分类专栏：数论/数学知识点 ACM中的数学问题合集 All

本文链接：https://blog.csdn.net/jk_chen_acmer/article/details/99889965

版权

All 同时被 3 个专栏收录

629 篇文章 6 订阅

订阅专栏

数论/数学知识点

48 篇文章 6 订阅

订阅专栏

ACM中的数学问题合集

44 篇文章 33 订阅

订阅专栏

Lagrange多项式插值

给出样例插值点，求出 $X$ 的函数值 $F (X)$ 。其中 $X\in R$

时间复杂度为 $O(k^2)$ ， $k$ 为插值点个数

连续函数：

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
type lagrange(vector<type> x,vector<type> y,type X)
{
	int n=x.size()-1;
	type ans=0;
	rep(k,0,n)
	{
		type temp=y[k];
		rep(j,0,n)if(j!=k)temp*=(X-x[j])/(x[k]-x[j]);
		ans+=temp;
	}
	return ans;
}
int main()
{
	vector<type> x;
	vector<type> y;
	rep(i,0,4)
        x.push_back(i),
        y.push_back(1.0*i*i+1.0*i/6);
	type X;
	while(cin>>X)cout<<lagrange(x,y,X)<<endl;
	return 0;
}
/*
f(x) = x*x + x/6
    0.5 -> 0.333333
    1.2 -> 1.64
*/

在取模的基础上同样适用，只需要将除法改成逆元即可。

例题为2018南京G题。（时间复杂度刚好，多了逆元会变成 $O(k^2logN)$ ，所以要先预处理逆元。跑了1200ms。离散函数）

#include <bits/stdc++.h>
#define rep(i,a,b) for(LL i=a;i<=b;i++)
using namespace std;
typedef long long LL;
typedef long long type;
const LL mod=1e9+7;
LL Pow(LL a,LL b,LL mod){
    LL res=1;
    while(b>0){
        if(b&1)res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}
LL inv[10][10];
type lagrange(vector<type> x,vector<type> y,type X)
{
	int n=x.size()-1;
	type ans=0;
	rep(i,0,n)
	{
		type temp=y[i];
		rep(j,0,n)if(j!=i)temp=(X-x[j])*inv[i][j]%mod*temp%mod;
		ans+=temp;
	}
	return (ans%mod+mod)%mod;
}
int main()
{
	vector<type> x={0,1,2,3,4};
	vector<type> y={0,1,5,15,35};
	int n=x.size()-1;
	rep(i,0,n)
        rep(j,0,n)
            inv[i][j]=Pow(x[i]-x[j],mod-2,mod );
	type X;
	int t;scanf("%d",&t);
	while(t--)scanf("%lld",&X),printf("%lld\n",lagrange(x,y,X));
	return 0;
}

使用Lagrange插值法求出多项式通项公式：（可能不准~）

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll gcd(ll a, ll b)
{
	return b ? gcd(b, a%b) : a;
}
class frac
{
	public:
	ll x,y;
	frac(){}frac(ll x,ll y):x(x),y(y){}
	bool operator < (const frac &b)const{return x*b.y<y*b.x;}
	bool operator > (const frac &b)const{return x*b.y>y*b.x;}
	bool operator ==(const frac &b)const{return x*b.y==y*b.x;}
	frac operator + (const frac &b)const{ll d=gcd(x*b.y+b.x*y,y*b.y);return frac((x*b.y+b.x*y)/d,(y*b.y)/d);}
	frac operator - (const frac &b)const{ll d=gcd(x*b.y-b.x*y,y*b.y);return frac((x*b.y-b.x*y)/d,(y*b.y)/d);}
	frac operator * (const frac &b)const{ll d=gcd(x*b.x,y*b.y);return frac((x*b.x)/d,(y*b.y)/d);}
	frac operator / (const frac &b)const{ll d=gcd(x*b.y,b.x*y);return frac((x*b.y)/d,(b.x*y)/d);}
	frac operator * (ll b)const{ll d=gcd(x*b,y);return frac((x*b)/d,(y)/d);}
	frac operator / (ll b)const{ll d=gcd(x,y*b);return frac((x)/d,(y*b)/d);}
	frac operator = (ll b){*this=frac(b,1);return *this;}
};
ostream &operator <<(ostream &out,const frac &a)
{
	if(a.y==1)out<<a.x;
	else out<<a.x<<"/"<<a.y;
	return out;
}
typedef frac type;
bool isZero(type x){
	return x.x==0;
}
class Poly{
public:
	vector<type>a={frac(0,1)};
	Poly(){}
	Poly(vector<type> b):a(b){}
	ll n(){
		return a.size()-1;
	}
	Poly operator = (type b){
		this->a.resize(1);
		this->a[0]=b;
		return *this;
	}
	Poly operator = (vector<type> b){
		this->a=b;
		return *this;
	}
	friend ostream &operator << (ostream &o,const Poly &f){
		for(int i=f.a.size()-1;~i;i--){
			if(!i)cout<<"("<<f.a[i]<<")";
			else cout<<"("<<f.a[i]<<")"<<"x^"<<i<<"+";
		}
		cout<<endl;
	}
	type coef(int i){
		if(i>=a.size() || i<0)return frac(0,1);
		return a[i];
	}
	type& operator [] (int i){
		if(i>=a.size() || i<0)cout<<" Warning: Index out of range\n";
		return a[i];
	}
	type operator () (type x){
		type ans;
		ans=0;
		for(int i=n();~i;i--)ans=ans*x+a[i];
		return ans;
	}
	Poly operator () (Poly x){
		Poly ans,t;
		for(int i=n();~i;i--){
			t=Poly((vector<type>){a[i]});
			ans=ans*x+t;
		}
		return ans;
	}
	Poly operator + (Poly &b){
		Poly c;
		c.a.resize(max(a.size(),b.a.size()));
		for(int i=0;i<c.a.size();i++)c.a[i]=coef(i)+b.coef(i);
		while(c.a.size()>1 && isZero(*(c.a.end()-1)))c.a.erase(c.a.end()-1);
		return c;
	}
	Poly operator - (Poly &b){
		Poly c;
		c.a.resize(max(a.size(),b.a.size()));
		for(int i=0;i<c.a.size();i++)c.a[i]=coef(i)-b.coef(i);
		while(c.a.size()>1 && isZero(*(c.a.end()-1)))c.a.erase(c.a.end()-1);
		return c;
	}
	Poly operator * (Poly &b){
		Poly c;
		c.a.resize(a.size()+b.a.size()-1);
		for(int i=0;i<c.a.size();i++)c.a[i]=0;
		for(int i=0;i<a.size();i++)
			for(int j=0;j<b.a.size();j++)
				c.a[i+j]=c.a[i+j]+a[i]*b.a[j];
		while(c.a.size()>1 && isZero(*(c.a.end()-1)))c.a.erase(c.a.end()-1);
		return c;
	}
};
Poly Lagrange(vector<type> x,vector<type> y){
	int n=x.size()-1;
	Poly ans;
	rep(k,0,n){
		Poly t,p;
		t=y[k];
		rep(j,0,n)if(j!=k){
			p=(vector<type>){frac(0,1)-x[j]/(x[k]-x[j]),frac(1,1)/(x[k]-x[j])};
			t=t*p;
		}
		ans=ans+t;
	}
	return ans;
}

int main(){
	vector<type> x={frac(0,1),frac(1,1),frac(2,1),frac(3,1),frac(4,1)};
	vector<type> y={frac(0,1),frac(1,1),frac(5,1),frac(15,1),frac(35,1)};
	Poly f=Lagrange(x,y);
	cout<<f;
	ll X;
	while(cin>>X)cout<<f(frac(X,1))<<endl;
	return 0;
}

Newton多项式插值

时间复杂度较为优秀，为 $O (k)$ 。连续函数：

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
class NewtonPoly{
	public:
	type f[105],d[105],x[105];
	ll n=0;
	void add(type X,type Y){
		x[n]=X,f[n]=Y;
		rep(i,1,n)f[n-i]=(f[n-i+1]-f[n-i])/(x[n]-x[n-i]);
		d[n++]=f[0];
	}
	type cal(type X){
		type ans=0,t=1;
		rep(i,0,n-1)ans+=d[i]*t,t*=X-x[i];
		return ans;
	}
}P;
int main(){
	P.add(0,0);
	P.add(1,1);
	P.add(2,5);
	P.add(3,15);
	P.add(4,35);
	type x;
	while(cin>>x)cout<<P.cal(x)<<endl;
	return 0;
}

南京G题AC代码（600ms，离散函数）：

#include <bits/stdc++.h>
#define rep(i,a,b) for(LL i=a;i<=b;i++)
using namespace std;
typedef long long LL;
typedef long long type;
const LL mod=1e9+7;
LL Pow(LL a,LL b,LL mod){
    LL res=1;
    while(b>0){
        if(b&1)res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}
class NewtonPoly{
	public:
	type f[105],d[105],x[105];
	LL n=0;
	void add(type X,type Y){
		x[n]=X,f[n]=Y;
		rep(i,1,n)f[n-i]=(f[n-i+1]-f[n-i])*Pow(x[n]-x[n-i],mod-2,mod)%mod;
		d[n++]=f[0];
	}
	type cal(type X){
		type ans=0,t=1;
		rep(i,0,n-1)ans=(ans+d[i]*t)%mod,t=t*(X-x[i])%mod;
		return (ans+mod)%mod;
	}
}P;
int main(){
	P.add(0,0);
	P.add(1,1);
	P.add(2,5);
	P.add(3,15);
	P.add(4,35);
	type x;
	int t;scanf("%d",&t);
	while(t--)scanf("%lld",&x),printf("%lld\n",P.cal(x));
	return 0;
}

高维插值

同样可以改成离散函数写法

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
typedef long double type;
type lagrange2(vector<type> x,vector<type> y,vector<vector<type> > z,type X,type Y){
	int M=x.size()-1,N=y.size()-1;
	type ans=0;
	rep(m,0,M)rep(n,0,N){
		type t=z[m][n];
		rep(i,0,M)if(i!=m)t*=(X-x[i])/(x[m]-x[i]);
		rep(i,0,N)if(i!=n)t*=(Y-y[i])/(y[n]-y[i]);
		ans+=t;
	}
	return ans;
}
int main(){
	vector<ll> x={0,1,2,3,4,5,6,7,8,9};
	vector<ll> y={0,1,2,3,4,5,6,7,8,9};
	vector<vector<ll> > z={
		{0,0,0,0,0,0,0,0,0,0},
		{0,1,2,3,4,5,6,7,8,9},
		{-2,2,6,10,14,18,22,26,30,34},
		{-10,0,10,20,30,40,50,60,70,80},
		{-30,-10,10,30,50,70,90,110,130,150},
		{-70,-35,0,35,70,105,140,175,210,245},
		{-140,-84,-28,28,84,140,196,252,308,364},
		{-252,-168,-84,0,84,168,252,336,420,504},
		{-420,-300,-180,-60,60,180,300,420,540,660},
		{-660,-495,-330,-165,0,165,330,495,660,825}
	};
	type X,Y;
	while(cin>>X>>Y)cout<<lagrange2(x,y,z,X,Y)<<endl;
	return 0;
}