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

/* mtst.c
Consistency tests for math functions.
Inverse function tests and Wronksians for random arguments.
With NTRIALS=10000, the following are typical results for
IEEE double precision arithmetic:

Consistency test of math functions.
Max and rms relative errors for 10000 random arguments.
x = cbrt( cube(x) ): max = 1.54E-016 rms = 6.50E-018
x = atan( tan(x) ): max = 2.83E-016 rms = 7.68E-017
x = sin( asin(x) ): max = 4.44E-016 rms = 7.01E-017
x = sqrt( square(x) ): max = 1.57E-016 rms = 2.70E-017
x = log( exp(x) ): max = 1.88E-014 rms = 2.02E-016
x = tanh( atanh(x) ): max = 4.00E-016 rms = 8.29E-017
x = asinh( sinh(x) ): max = 2.04E-016 rms = 1.56E-017
x = acosh( cosh(x) ): max = 2.67E-014 rms = 2.95E-016
x = pow( pow(x,a),1/a ): max = 3.38E-012 rms = 3.39E-014
Legendre ellpk, ellpe: max = 1.55E-011 rms = 1.74E-013
x = ndtri( ndtr(x) ): max = 4.09E-014 rms = 6.40E-016
Absolute error criterion (but relative if >1):
lgam(x) = log(gamma(x)): max = 4.49E-016 rms = 8.80E-017
x = cos( acos(x) ): max = 4.44E-016 rms = 1.00E-016
Absolute error and only 2000 trials:
Wronksian of Yn, Jn: max = 7.12E-011 rms = 1.59E-012
Wronksian of Kn, Iv: max = 7.07E-012 rms = 5.25E-013

*/

/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include "mconf.h"

#define NTRIALS 1000
#define WTRIALS (NTRIALS/5)
#define STRTST 0

double ndtr(), ndtri(), ellpe(), ellpk(), gamma(), lgam();
double fabs(), sqrt();
double cbrt(), exp(), log(), tan(), atan();
double sin(), asin(), cos(), acos(), pow();
double tanh(), atanh(), sinh(), asinh(), cosh(), acosh();
double jn(), yn(), iv(), kn();

/* Provide inverses for square root and cube root: */
double square(x)
double x;
{
return( x * x );
}

double cube(x)
double x;
{
return( x * x * x );
}

/* lookup table for each function */
struct fundef
{
char *nam1; /* the function */
double (*name )();
char *nam2; /* its inverse */
double (*inv )();
int nargs; /* number of function arguments */
int tstyp; /* type code of the function */
long ctrl; /* relative error flag */
double arg1w; /* width of domain for 1st arg */
double arg1l; /* lower bound domain 1st arg */
long arg1f; /* flags, e.g. integer arg */
double arg2w; /* same info for args 2, 3, 4 */
double arg2l;
long arg2f;
/*
double arg3w;
double arg3l;
long arg3f;
double arg4w;
double arg4l;
long arg4f;
*/
};


/* fundef.ctrl bits: */
#define RELERR 1
#define EXPSCAL 4

/* fundef.tstyp test types: */
#define POWER 1
#define ELLIP 2
#define GAMMA 3
#define WRONK1 4
#define WRONK2 5
#define WRONK3 6

/* fundef.argNf argument flag bits: */
#define INT 2

extern double MINLOG;
extern double MAXLOG;
extern double PI;
extern double PIO2;
/*
define MINLOG -170.0
define MAXLOG +170.0
define PI 3.14159265358979323846
define PIO2 1.570796326794896619
*/

#define NTESTS 15
struct fundef defs[NTESTS] = {
{" cube", cube, " cbrt", cbrt, 1, 0, 1, 2000.0, -1000.0, 0,
0.0, 0.0, 0},
{" tan", tan, " atan", atan, 1, 0, 1, 0.0, 0.0, 0,
0.0, 0.0, 0},
{" asin", asin, " sin", sin, 1, 0, 1, 2.0, -1.0, 0,
0.0, 0.0, 0},
{"square", square, " sqrt", sqrt, 1, 0, 1, 170.0, -85.0, EXPSCAL,
0.0, 0.0, 0},
{" exp", exp, " log", log, 1, 0, 1, 340.0, -170.0, 0,
0.0, 0.0, 0},
{" atanh", atanh, " tanh", tanh, 1, 0, 1, 2.0, -1.0, 0,
0.0, 0.0, 0},
{" sinh", sinh, " asinh", asinh, 1, 0, 1, 340.0, 0.0, 0,
0.0, 0.0, 0},
{" cosh", cosh, " acosh", acosh, 1, 0, 1, 340.0, 0.0, 0,
0.0, 0.0, 0},
{"pow", pow, "pow", pow, 2, POWER, 1, 25.0, 0.0, 0,
50.0, -25.0, 0},
{" ellpe", ellpe, " ellpk", ellpk, 1, ELLIP, 1, 1.0, 0.0, 0,
0.0, 0.0, 0},
{" ndtr", ndtr, " ndtri", ndtri, 1, 0, 1, 10.0, -10.0, 0,
0.0, 0.0, 0},
/* Absolute error criterion starts with tstyp = GAMMA: */
{"gamma", gamma, "lgam", lgam, 1, GAMMA, 0, 11.0, 1.0, 0,
0.0, 0.0, 0},
{" acos", acos, " cos", cos, 1, 0, 0, 2.0, -1.0, 0,
0.0, 0.0, 0},
/* Smaller number of trials starts with tstyp = WRONK1: */
{" Jn", jn, " Yn", yn, 2, WRONK1, 0, 30.0, 0.1, 0,
40.0, -20.0, INT},
{" Iv", iv, " Kn", kn, 2, WRONK2, 0, 30.0, 0.1, 0,
20.0, 0.0, INT},
};

static char *headrs[] = {
"x = %s( %s(x) ): ",
"x = %s( %s(x,a),1/a ): ", /* power */
"Legendre %s, %s: ", /* ellip */
"%s(x) = log(%s(x)): ", /* gamma */
"Wronksian of %s, %s: ",
"Wronksian of %s, %s: ",
"Wronksian of %s, %s: "
};

static double y1 = 0.0;
static double y2 = 0.0;
static double y3 = 0.0;
static double y4 = 0.0;
static double a = 0.0;
static double b = 0.0;
static double c = 0.0;
static double x = 0.0;
static double y = 0.0;
static double z = 0.0;
static double e = 0.0;
static double max = 0.0;
static double rmsa = 0.0;
static double rms = 0.0;
static double ave = 0.0;


main()
{
double (*fun )();
double (*ifun )();
struct fundef *d;
int i, j, k, itst;
int m, ntr;

ntr = NTRIALS;
printf( "Consistency test of math functions.\n" );
printf( "Max and rms relative errors for %d random arguments.\n",
ntr );

/* Initialize machine dependent parameters: */
defs[1].arg1w = PI;
defs[1].arg1l = -PI/2.0;
defs[3].arg1w = MAXLOG;
defs[3].arg1l = -MAXLOG/2.0;
defs[4].arg1w = 2.0*MAXLOG;
defs[4].arg1l = -MAXLOG;
defs[6].arg1w = 2.0*MAXLOG;
defs[6].arg1l = -MAXLOG;
defs[7].arg1w = MAXLOG;
defs[7].arg1l = 0.0;


/* Outer loop, on the test number: */

for( itst=STRTST; itst {
d = &defs[itst];
m = 0;
max = 0.0;
rmsa = 0.0;
ave = 0.0;
fun = d->name;
ifun = d->inv;

/* Absolute error criterion starts with gamma function
* (put all such at end of table)
*/
if( d->tstyp == GAMMA )
printf( "Absolute error criterion (but relative if >1):\n" );

/* Smaller number of trials for Wronksians
* (put them at end of list)
*/
if( d->tstyp == WRONK1 )
{
ntr = WTRIALS;
printf( "Absolute error and only %d trials:\n", ntr );
}

printf( headrs[d->tstyp], d->nam2, d->nam1 );

for( i=0; i {
m++;

/* make random number(s) in desired range(s) */
switch( d->nargs )
{

default:
goto illegn;

case 2:
drand( &a );
a = d->arg2w * ( a - 1.0 ) + d->arg2l;
if( d->arg2f & EXPSCAL )
{
a = exp(a);
drand( &y2 );
a -= 1.0e-13 * a * y2;
}
if( d->arg2f & INT )
{
k = a + 0.25;
a = k;
}

case 1:
drand( &x );
x = d->arg1w * ( x - 1.0 ) + d->arg1l;
if( d->arg1f & EXPSCAL )
{
x = exp(x);
drand( &a );
x += 1.0e-13 * x * a;
}
}


/* compute function under test */
switch( d->nargs )
{
case 1:
switch( d->tstyp )
{
case ELLIP:
y1 = ( *(fun) )(x);
y2 = ( *(fun) )(1.0-x);
y3 = ( *(ifun) )(x);
y4 = ( *(ifun) )(1.0-x);
break;

case GAMMA:
y = lgam(x);
x = log( gamma(x) );
break;

default:
z = ( *(fun) )(x);
y = ( *(ifun) )(z);
}
break;

case 2:
if( d->arg2f & INT )
{
switch( d->tstyp )
{
case WRONK1:
y1 = (*fun)( k, x ); /* jn */
y2 = (*fun)( k+1, x );
y3 = (*ifun)( k, x ); /* yn */
y4 = (*ifun)( k+1, x );
break;

case WRONK2:
y1 = (*fun)( a, x ); /* iv */
y2 = (*fun)( a+1.0, x );
y3 = (*ifun)( k, x ); /* kn */
y4 = (*ifun)( k+1, x );
break;

default:
z = (*fun)( k, x );
y = (*ifun)( k, z );
}
}
else
{
if( d->tstyp == POWER )
{
z = (*fun)( x, a );
y = (*ifun)( z, 1.0/a );
}
else
{
z = (*fun)( a, x );
y = (*ifun)( a, z );
}
}
break;


default:
illegn:
printf( "Illegal nargs= %d", d->nargs );
exit(1);
}

switch( d->tstyp )
{
case WRONK1:
e = (y2*y3 - y1*y4) - 2.0/(PI*x); /* Jn, Yn */
break;

case WRONK2:
e = (y2*y3 + y1*y4) - 1.0/x; /* In, Kn */
break;

case ELLIP:
e = (y1-y3)*y4 + y3*y2 - PIO2;
break;

default:
e = y - x;
break;
}

if( d->ctrl & RELERR )
e /= x;
else
{
if( fabs(x) > 1.0 )
e /= x;
}

ave += e;
/* absolute value of error */
if( e < 0 )
e = -e;

/* peak detect the error */
if( e > max )
{
max = e;

if( e > 1.0e-10 )
{
printf("x %.6E z %.6E y %.6E max %.4E\n",
x, z, y, max);
if( d->tstyp >= WRONK1 )
{
printf( "y1 %.4E y2 %.4E y3 %.4E y4 %.4E k %d x %.4E\n",
y1, y2, y3, y4, k, x );
}
}

/*
printf("%.8E %.8E %.4E %6ld \n", x, y, max, n);
printf("%d %.8E %.8E %.4E %6ld \n", k, x, y, max, n);
printf("%.6E %.6E %.6E %.4E %6ld \n", a, x, y, max, n);
printf("%.6E %.6E %.6E %.6E %.4E %6ld \n", a, b, x, y, max, n);
printf("%.4E %.4E %.4E %.4E %.4E %.4E %6ld \n",
a, b, c, x, y, max, n);
*/
}

/* accumulate rms error */
e *= 1.0e16; /* adjust range */
rmsa += e * e; /* accumulate the square of the error */
endlup:
;
}

/* report after NTRIALS trials */
rms = 1.0e-16 * sqrt( rmsa/m );
printf(" max = %.2E rms = %.2E\n", max, rms );
} /* loop on itst */

}