/* remesf.c */

#include "remes.h"


/* The flag bits for type of approximation: */
/* PXSQ | XPX | X2PX | SQL | SQH | PADE | CW | ZER | SPECIAL | PXCU */
int config =   ZER ;

/* Insert function name and formulas for printout */
char funnam[] = {
 "log(1+x) "
};

char znam[] = { "x" };
double sqrt(), log(), exp(), pow(), sin(), cos();
double j0(), y0(), atan(), floor();
double j1(), y1();
double exp10(), log10();
double erfc();
double lx();

/* This subroutine computes the rational form P(x)/Q(x) */
/* using coefficients from the solution vector param[]. */

double approx(x)
double x;
{
double gx, z, yn, yd, q;
double gofx(), speci();
int i;

gx = gofx(x);
if( config & PXCU )
	z = gx * gx * gx;
else if( config & PXSQ )
	z = gx * gx;
else
	z = gx;

/* Highest order numerator coefficient */
yn = param[n];

/* Work backwards toward the constant term. */
for( i=n-1; i>=0; i-- )
	yn = z * yn  +  param[i];
	
if( d > 0 )
	{
/* Highest degree coefficient = 1.0 */
	yd = z + param[n+d];

	for( i=n+d-1; i>n; i-- )
		yd = z * yd  +  param[i];
	}
else
/* There is no denominator. */
	yd = 1.0;

if( config & XPX )
	yn = yn * gx;
if( config & X2PX )
	yn = yn * gx * gx;
if( config & PADE )
	{ /* 2P/(Q-P) */
	yd = yd - yn;
	yn = 2.0 * yn;
	}
qyaprx = yn/yd;
if( config & CW )
	qyaprx = gx + qyaprx * gx * gx;
if( config & SPECIAL )
	qyaprx = speci( qyaprx, gx );
return( qyaprx );
}



/* Subroutine to compute approximation error at x */

double geterr(x)
double x;
{
double e, f;
double fabs(), approx(), func();

f = func(x);
e = approx(x) - f;
if( relerr )
	{
	if( f != 0.0 )
		e /= f;
	}
if( e < 0.0 )
	{
	esign = -1;
	e = -e;
	}
else
	esign = 1;

return(e);
}



/* Subroutine for special argument transformations */

double gofx(x)
double x;
{

return(x);
}

/* Routine for special modifications of the approximating form.
 * Example already provided by the CW flag: 
 *	y(1+dev) = gx + gx^2 R(gx)
 * would change y to
 *	R(gx) = (y - gx)/(gx*gx)
 * This function is called from remese.c.
 *
 * An inverse routine called speci() must also be supplied.
 * This finds y from R and gx (see below).
 */
extern double PI, PIO4, THPIO4;
#define LS2PI 0.91893853320467274178
#define SQTPI 2.50662827463100050242
#define JZ1   5.783185962946784521176
#define JZ2  30.471262343662086399078
#define JZ3  74.887006790695183444889
#define YZ1  0.43221455686510834878
#define YZ2 22.401876406482861405
#define YZ3 64.130620282338755553298
#define JO1 14.6819706421238932572
#define YO1  4.66539330185668857532


extern double TWOOPI;

double special( y, gx )
double y, gx;
{
double a, b;
/* exponential, y = exp(x) = 1 + x + .5x^2 + x^2 R(x) */
/*
if( gx == 0.0 )
	return(1.0);
b = gx * gx;
a = (y - 1.0 - gx - .5*b)/b;
*/

/* y = exp2(x) = 1 + x R(x) */
/*
a = y - 1.0;
*/

/* y = cos(x) = 1 - .5 x^2 + x^2 x^2 P(x^2) */
/*
b = gx * gx;
a = (y - 1.0 + 0.5*b)/b;
*/

/* logarithm, y = log(1+x) = x - .5x^2 + x^2 * (xP(x))
 * configuration is  ZER | XPX | SPECIAL
 */
/*
if( gx == 0.0 )
	return(0.0);
b = gx * gx;
a = (y - gx + 0.5 * b)/b;
*/

/* y = log gamma(x) = q(x) + 1/x P(1/x^2)
 * configufation is SQH | ZER | XPX | PXSQ | SPECIAL
 */
/*
b = 1.0/gx;
b = ( b - 0.5 ) * log(b) - b + LS2PI;
a = y - b;
*/

/* y = gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) */
/*
b = 1.0/gx;
a = SQTPI * pow( b, b-0.5 ) * exp(-b);
a = (y - a)/a;
*/

/* y = j0(x) = (x^2 - JZ1)(x^2-JZ2)(x^2-JZ3)P(x^2) */
/*
b=gx*gx;
a = (b-JZ1);
a = y/a;
*/

/* y = y0(x) = TWOOPI * log(x) * j0(x) + (x^2-YZ1)*P(x^2) */
/*
b=gx*gx;
a = y - TWOOPI * log(gx) * j0(gx);
a /= (b-YZ1);
*/

/* y = j1(x) = (x^2 - JO1)P(x^2) */
/*
b=gx*gx;
a = (b-JO1);
a = y/a;
*/
/* y = y1(x) = TWOOPI * (log(x) * j1(x) - 1/x) + (x^2-YZ1)*P(x^2) */
/*
b = gx*gx;
a = y - TWOOPI * ( j1(gx) * log(gx)  -  1.0/gx );
a /= (b-YO1);
*/

/* y = exp2(x) = 1 + x P(x) */
a = y - 1.0;

/* y = erfc(x) = exp(-x^2) P(x) */
a = y * exp( gx * gx );

/* acosh() */
/*
if( gx == 0.0 )
	return(0.0);
a = y/(2.0*sqrt(gx));
*/
return( a );
}



double speci( R, gx )
double R, gx;
{
double y, b;

/* exponential, y = exp(x) = 1 + x + .5x^2 + x^2 R(x) */
/*
b = gx * gx;
y = 1.0 + gx + .5 * b;
y = y + b * R;
*/

/* y = exp2(x) = 1 + x R(x) */
/*
y = R + 1.0;
*/
/* y = cos(x) = 1 - .5 x^2 +  x^2 x^2 R(x^2) */
/*
b = gx * gx;
y = 1.0 - 0.5*b + b * R;
*/

/* log(1+x) = x - .5x^2 + x^2 xR(x) */
/*
b = gx * gx;
y = gx  - 0.5 * b  +  b * R;
*/

/* y = log gamma(x) = q(x) + 1/x P(1/x^2) */
/*
b = 1.0/gx;
b = ( b - 0.5 ) * log(b) - b + LS2PI;
y = b + R;
*/

/* y = gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) */
/*
b = 1.0/gx;
b = SQTPI * pow( b, b-0.5 ) * exp(-b);
y = b + b * R;
*/
/* y = j0(x) = (x^2 - JZ1)(x^2-JZ2)(x^2-JZ3)P(x^2) */
/*
b=gx*gx;
y = (b-JZ1)*R;
*/

/* y = y0(x) = TWOOPI * log(x) * j0(x) + (x^2-YZ1)*P(x^2) */
/*
b=gx*gx;
y = TWOOPI * log(gx) * j0(gx) + (b-YZ1)*R;
*/

/* y = j1(x) = (x^2 - JO1)P(x^2) */
/*
b=gx*gx;
y = (b-JO1)*R;
*/

/* y = y1(x) = TWOOPI * (log(x) * j1(x) - 1/x) + (x^2-YZ1)*P(x^2) */
/*
b = gx*gx;
y =  TWOOPI * ( j1(gx) * log(gx)  -  1.0/gx ) + (b-YO1)*R;
*/

/* y = exp2(x) = 1 + x P(x) */
/*y = 1.0 + R;*/

/* y = erfc(x) = exp(-x^2) P(x) */
y = exp( -gx * gx ) * R;

/*y = 2.0 * sqrt(gx) * R;*/
return( y );
}

/* Put here an accurate routine */
/* for the function to be approximated. */

static int fflg = 0;
static double ff = 0.0;

double func(x)
double x;
{
double xx, y, t, u, s, c;


qx = x;

if( fflg == 0 )
	{
	fflg = 1;
	ff = 10.0 * log10(2.0);
	}


if( x == 0.0 )
	{
/*	y = ff;*/
	y = 0.0;
	goto done;
	}

/* Bessel, phase */
/*
xx = 1.0/x;
y = j0(xx);
t = y0(xx);
y = atan(t/y);
t = xx - PIO4;
t = t - PI * floor(t/PI + 0.5);
y -= t;
if( y > 0.5*PI )
	y -= PI;
if( y < -0.5*PI )
	y += PI;
*/
/* Bessel, modulus */
/*
xx = 1.0/(x*x);
y = j1(xx);
t = y1(xx);
y = sqrt( y*y + t*t );
*/

/* bes1, phase */
/*
xx = 1.0/x;
y = j1(xx);
t = y1(xx);
y = atan(t/y);
t = xx - THPIO4;
t = t - PI * floor(t/PI + 0.5);
y -= t;
if( y > 0.5*PI )
	y -= PI;
if( y < -0.5*PI )
	y += PI;
*/

y = log(1.0+x);

done:
qy = y;
return( y );
}