Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : REMESF.C

/* 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 );
}