From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: Totally rework the math library, this time based on the MacOs X math library (which is itself based on the math lib from FreeBSD). -Erik --- libm/double/mtst.c | 464 ----------------------------------------------------- 1 file changed, 464 deletions(-) delete mode 100644 libm/double/mtst.c (limited to 'libm/double/mtst.c') diff --git a/libm/double/mtst.c b/libm/double/mtst.c deleted file mode 100644 index 2559d2340..000000000 --- a/libm/double/mtst.c +++ /dev/null @@ -1,464 +0,0 @@ -/* mtst.c - Consistency tests for math functions. - To get strict rounding rules on a 386 or 68000 computer, - define SETPREC to 1. - - 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 = 0.00E+00 rms = 0.00E+00 -x = atan( tan(x) ): max = 2.21E-16 rms = 3.27E-17 -x = sin( asin(x) ): max = 2.13E-16 rms = 2.95E-17 -x = sqrt( square(x) ): max = 0.00E+00 rms = 0.00E+00 -x = log( exp(x) ): max = 1.11E-16 A rms = 4.35E-18 A -x = tanh( atanh(x) ): max = 2.22E-16 rms = 2.43E-17 -x = asinh( sinh(x) ): max = 2.05E-16 rms = 3.49E-18 -x = acosh( cosh(x) ): max = 1.43E-15 A rms = 1.54E-17 A -x = log10( exp10(x) ): max = 5.55E-17 A rms = 1.27E-18 A -x = pow( pow(x,a),1/a ): max = 7.60E-14 rms = 1.05E-15 -x = cos( acos(x) ): max = 2.22E-16 A rms = 6.90E-17 A -*/ - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier -*/ - - -#include -#include -#include - -#ifndef NTRIALS -#define NTRIALS 10000 -#endif - -#define SETPREC 1 -#define STRTST 0 - -#define WTRIALS (NTRIALS/5) - -#ifdef ANSIPROT -extern double fabs ( double ); -extern double sqrt ( double ); -extern double cbrt ( double ); -extern double exp ( double ); -extern double log ( double ); -extern double exp10 ( double ); -extern double log10 ( double ); -extern double tan ( double ); -extern double atan ( double ); -extern double sin ( double ); -extern double asin ( double ); -extern double cos ( double ); -extern double acos ( double ); -extern double pow ( double, double ); -extern double tanh ( double ); -extern double atanh ( double ); -extern double sinh ( double ); -extern double asinh ( double x ); -extern double cosh ( double ); -extern double acosh ( double ); -extern double gamma ( double ); -extern double lgam ( double ); -#else -double fabs(), sqrt(), cbrt(), exp(), log(); -double exp10(), log10(), tan(), atan(); -double sin(), asin(), cos(), acos(), pow(); -double tanh(), atanh(), sinh(), asinh(), cosh(), acosh(); -double gamma(), lgam(); -#endif - -/* C9X spells lgam lgamma. */ -#define GLIBC2 0 -#if GLIBC2 -double lgamma (double); -#endif - -#if SETPREC -int dprec(); -#endif - -int drand(); -/* void exit(); */ -/* int printf(); */ - - -/* 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 - -/* 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 -#define EXPSCAL 4 - -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 12 -struct fundef defs[NTESTS] = { -{" cube", cube, " cbrt", cbrt, 1, 0, 1, 2002.0, -1001.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, 0, 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, 0, 340.0, 0.0, 0, -0.0, 0.0, 0}, -{" exp10", exp10, " log10", log10, 1, 0, 0, 340.0, -170.0, 0, -0.0, 0.0, 0}, -{"pow", pow, "pow", pow, 2, POWER, 1, 21.0, 0.0, 0, -42.0, -21.0, 0}, -{" acos", acos, " cos", cos, 1, 0, 0, 2.0, -1.0, 0, -0.0, 0.0, 0}, -#if GLIBC2 -{ "gamma", gamma, "lgamma", lgamma, 1, GAMMA, 0, 34.0, 0.0, 0, -0.0, 0.0, 0}, -#else -{ "gamma", gamma, "lgam", lgam, 1, GAMMA, 0, 34.0, 0.0, 0, -0.0, 0.0, 0}, -#endif -}; - -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 yy1 = 0.0; -static double y2 = 0.0; -static double y3 = 0.0; -static double y4 = 0.0; -static double a = 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; - - -int main() -{ -double (*fun )(); -double (*ifun )(); -struct fundef *d; -int i, k, itst; -int m, ntr; - -#if SETPREC -dprec(); /* set coprocessor precision */ -#endif -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; -/* Microsoft C has trouble with denormal numbers. */ -#if 0 -defs[3].arg1w = MAXLOG; -defs[3].arg1l = -MAXLOG/2.0; -defs[4].arg1w = 2*MAXLOG; -defs[4].arg1l = -MAXLOG; -#endif -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; itstname; -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; inargs ) -{ - -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: - yy1 = ( *(fun) )(x); - y2 = ( *(fun) )(1.0-x); - y3 = ( *(ifun) )(x); - y4 = ( *(ifun) )(1.0-x); - break; - -#if 1 - case GAMMA: -#if GLIBC2 - y = lgamma(x); -#else - y = lgam(x); -#endif - x = log( gamma(x) ); - break; -#endif - default: - z = ( *(fun) )(x); - y = ( *(ifun) )(z); - } - break; - - case 2: - if( d->arg2f & INT ) - { - switch( d->tstyp ) - { - case WRONK1: - yy1 = (*fun)( k, x ); /* jn */ - y2 = (*fun)( k+1, x ); - y3 = (*ifun)( k, x ); /* yn */ - y4 = (*ifun)( k+1, x ); - break; - - case WRONK2: - yy1 = (*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 - yy1*y4) - 2.0/(PI*x); /* Jn, Yn */ - break; - - case WRONK2: - e = (y2*y3 + yy1*y4) - 1.0/x; /* In, Kn */ - break; - - case ELLIP: - e = (yy1-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 == POWER ) - { - printf( "a %.6E\n", a ); - } - if( d->tstyp >= WRONK1 ) - { - printf( "yy1 %.4E y2 %.4E y3 %.4E y4 %.4E k %d x %.4E\n", - yy1, 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 */ -} - -/* report after NTRIALS trials */ -rms = 1.0e-16 * sqrt( rmsa/m ); -if(d->ctrl & RELERR) - printf(" max = %.2E rms = %.2E\n", max, rms ); -else - printf(" max = %.2E A rms = %.2E A\n", max, rms ); -} /* loop on itst */ - -exit(0); -} -- cgit v1.2.3