Lagrange多项式插值
给出样例插值点,求出 X X X的函数值 F ( X ) F(X) F(X)。其中 X ∈ R X\in R X∈R
时间复杂度为 O ( k 2 ) O(k^2) O(k2), k k 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 2 l o g N ) O(k^2logN) O(k2logN),所以要先预处理逆元。跑了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 ) O(k) 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;
}