text: struct Base_interp
  {
  	Int n, mm, jsav, cor, dj;
  	const Doub *xx, *yy;
  	Base_interp(VecDoub_I &x, const Doub *y, Int m)
  		: n(x.size()), mm(m), jsav(0), cor(0), xx(&x[0]), yy(y) {
  		dj = MIN(1,(int)pow((Doub)n,0.25));
  	}
  
  	Doub interp(Doub x) {
  		Int jlo = cor ? hunt(x) : locate(x);
  		return rawinterp(jlo,x);
  	}
  
  	Int locate(const Doub x);
  	Int hunt(const Doub x);
  	
  	Doub virtual rawinterp(Int jlo, Doub x) = 0;
  
  };
  Int Base_interp::locate(const Doub x)
  {
  	Int ju,jm,jl;
  	if (n < 2 || mm < 2 || mm > n) throw("locate size error");
  	Bool ascnd=(xx[n-1] >= xx[0]);
  	jl=0;
  	ju=n-1;
  	while (ju-jl > 1) {
  		jm = (ju+jl) >> 1;
  		if (x >= xx[jm] == ascnd)
  			jl=jm;
  		else
  			ju=jm;
  	}
  	cor = abs(jl-jsav) > dj ? 0 : 1;
  	jsav = jl;
  	return MAX(0,MIN(n-mm,jl-((mm-2)>>1)));
  }
  Int Base_interp::hunt(const Doub x)
  {
  	Int jl=jsav, jm, ju, inc=1;
  	if (n < 2 || mm < 2 || mm > n) throw("hunt size error");
  	Bool ascnd=(xx[n-1] >= xx[0]);
  	if (jl < 0 || jl > n-1) {
  		jl=0;
  		ju=n-1;
  	} else {
  		if (x >= xx[jl] == ascnd) {
  			for (;;) {
  				ju = jl + inc;
  				if (ju >= n-1) { ju = n-1; break;}
  				else if (x < xx[ju] == ascnd) break;
  				else {
  					jl = ju;
  					inc += inc;
  				}
  			}
  		} else {
  			ju = jl;
  			for (;;) {
  				jl = jl - inc;
  				if (jl <= 0) { jl = 0; break;}
  				else if (x >= xx[jl] == ascnd) break;
  				else {
  					ju = jl;
  					inc += inc;
  				}
  			}
  		}
  	}
  	while (ju-jl > 1) {
  		jm = (ju+jl) >> 1;
  		if (x >= xx[jm] == ascnd)
  			jl=jm;
  		else
  			ju=jm;
  	}
  	cor = abs(jl-jsav) > dj ? 0 : 1;
  	jsav = jl;
  	return MAX(0,MIN(n-mm,jl-((mm-2)>>1)));
  }
  struct Poly_interp : Base_interp
  {
  	Doub dy;
  	Poly_interp(VecDoub_I &xv, VecDoub_I &yv, Int m)
  		: Base_interp(xv,&yv[0],m), dy(0.) {}
  	Doub rawinterp(Int jl, Doub x);
  };
  
  Doub Poly_interp::rawinterp(Int jl, Doub x)
  {
  	Int i,m,ns=0;
  	Doub y,den,dif,dift,ho,hp,w;
  	const Doub *xa = &xx[jl], *ya = &yy[jl];
  	VecDoub c(mm),d(mm);
  	dif=abs(x-xa[0]);
  	for (i=0;i<mm;i++) {
  		if ((dift=abs(x-xa[i])) < dif) {
  			ns=i;
  			dif=dift;
  		}
  		c[i]=ya[i];
  		d[i]=ya[i];
  	}
  	y=ya[ns--];
  	for (m=1;m<mm;m++) {
  		for (i=0;i<mm-m;i++) {
  			ho=xa[i]-x;
  			hp=xa[i+m]-x;
  			w=c[i+1]-d[i];
  			if ((den=ho-hp) == 0.0) throw("Poly_interp error");
  			den=w/den;
  			d[i]=hp*den;
  			c[i]=ho*den;
  		}
  		y += (dy=(2*(ns+1) < (mm-m) ? c[ns+1] : d[ns--]));
  	}
  	return y;
  }
  struct Rat_interp : Base_interp
  {
  	Doub dy;
  	Rat_interp(VecDoub_I &xv, VecDoub_I &yv, Int m)
  		: Base_interp(xv,&yv[0],m), dy(0.) {}
  	Doub rawinterp(Int jl, Doub x);
  };
  
  Doub Rat_interp::rawinterp(Int jl, Doub x)
  {
  	const Doub TINY=1.0e-99;
  	Int m,i,ns=0;
  	Doub y,w,t,hh,h,dd;
  	const Doub *xa = &xx[jl], *ya = &yy[jl];
  	VecDoub c(mm),d(mm);
  	hh=abs(x-xa[0]);
  	for (i=0;i<mm;i++) {
  		h=abs(x-xa[i]);
  		if (h == 0.0) {
  			dy=0.0;
  			return ya[i];
  		} else if (h < hh) {
  			ns=i;
  			hh=h;
  		}
  		c[i]=ya[i];
  		d[i]=ya[i]+TINY;
  	}
  	y=ya[ns--];
  	for (m=1;m<mm;m++) {
  		for (i=0;i<mm-m;i++) {
  			w=c[i+1]-d[i];
  			h=xa[i+m]-x;
  			t=(xa[i]-x)*d[i]/h;
  			dd=t-c[i+1];
  			if (dd == 0.0) throw("Error in routine ratint");
  			dd=w/dd;
  			d[i]=c[i+1]*dd;
  			c[i]=t*dd;
  		}
  		y += (dy=(2*(ns+1) < (mm-m) ? c[ns+1] : d[ns--]));
  	}
  	return y;
  }
  struct Spline_interp : Base_interp
  {
  	VecDoub y2;
  	
  	Spline_interp(VecDoub_I &xv, VecDoub_I &yv, Doub yp1=1.e99, Doub ypn=1.e99)
  	: Base_interp(xv,&yv[0],2), y2(xv.size())
  	{sety2(&xv[0],&yv[0],yp1,ypn);}
  
  	Spline_interp(VecDoub_I &xv, const Doub *yv, Doub yp1=1.e99, Doub ypn=1.e99)
  	: Base_interp(xv,yv,2), y2(xv.size())
  	{sety2(&xv[0],yv,yp1,ypn);}
  
  	void sety2(const Doub *xv, const Doub *yv, Doub yp1, Doub ypn);
  	Doub rawinterp(Int jl, Doub xv);
  };
  void Spline_interp::sety2(const Doub *xv, const Doub *yv, Doub yp1, Doub ypn)
  {
  	Int i,k;
  	Doub p,qn,sig,un;
  	Int n=y2.size();
  	VecDoub u(n-1);
  	if (yp1 > 0.99e99)
  		y2[0]=u[0]=0.0;
  	else {
  		y2[0] = -0.5;
  		u[0]=(3.0/(xv[1]-xv[0]))*((yv[1]-yv[0])/(xv[1]-xv[0])-yp1);
  	}
  	for (i=1;i<n-1;i++) {
  		sig=(xv[i]-xv[i-1])/(xv[i+1]-xv[i-1]);
  		p=sig*y2[i-1]+2.0;
  		y2[i]=(sig-1.0)/p;
  		u[i]=(yv[i+1]-yv[i])/(xv[i+1]-xv[i]) - (yv[i]-yv[i-1])/(xv[i]-xv[i-1]);
  		u[i]=(6.0*u[i]/(xv[i+1]-xv[i-1])-sig*u[i-1])/p;
  	}
  	if (ypn > 0.99e99)
  		qn=un=0.0;
  	else {
  		qn=0.5;
  		un=(3.0/(xv[n-1]-xv[n-2]))*(ypn-(yv[n-1]-yv[n-2])/(xv[n-1]-xv[n-2]));
  	}
  	y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
  	for (k=n-2;k>=0;k--)
  		y2[k]=y2[k]*y2[k+1]+u[k];
  }
  Doub Spline_interp::rawinterp(Int jl, Doub x)
  {
  	Int klo=jl,khi=jl+1;
  	Doub y,h,b,a;
  	h=xx[khi]-xx[klo];
  	if (h == 0.0) throw("Bad input to routine splint");
  	a=(xx[khi]-x)/h;
  	b=(x-xx[klo])/h;
  	y=a*yy[klo]+b*yy[khi]+((a*a*a-a)*y2[klo]
  		+(b*b*b-b)*y2[khi])*(h*h)/6.0;
  	return y;
  }
  struct BaryRat_interp : Base_interp
  {
  	VecDoub w;
  	Int d;
  	BaryRat_interp(VecDoub_I &xv, VecDoub_I &yv, Int dd);
  	Doub rawinterp(Int jl, Doub x);
  	Doub interp(Doub x);
  };
  
  BaryRat_interp::BaryRat_interp(VecDoub_I &xv, VecDoub_I &yv, Int dd)
  		: Base_interp(xv,&yv[0],xv.size()), w(n), d(dd)
  {
  	if (n<=d) throw("d too large for number of points in BaryRat_interp");
  	for (Int k=0;k<n;k++) {
  		Int imin=MAX(k-d,0);
  		Int imax = k >= n-d ? n-d-1 : k;
  		Doub temp = imin & 1 ? -1.0 : 1.0;
  		Doub sum=0.0;
  		for (Int i=imin;i<=imax;i++) {
  			Int jmax=MIN(i+d,n-1);
  			Doub term=1.0;
  			for (Int j=i;j<=jmax;j++) {
  				if (j==k) continue;
  				term *= (xx[k]-xx[j]);
  			}
  			term=temp/term;
  			temp=-temp;
  			sum += term;
  		}
  		w[k]=sum;
  	}
  }
  Doub BaryRat_interp::rawinterp(Int jl, Doub x)
  {
  	Doub num=0,den=0;
  	for (Int i=0;i<n;i++) {
  		Doub h=x-xx[i];
  		if (h == 0.0) {
  			return yy[i];
  		} else {
  			Doub temp=w[i]/h;
  			num += temp*yy[i];
  			den += temp;
  		}
  	}
  	return num/den;
  }
  Doub BaryRat_interp::interp(Doub x) {
  	return rawinterp(1,x);
  }
  