summaryrefslogtreecommitdiff
path: root/libm/double/mtst.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/mtst.c')
-rw-r--r--libm/double/mtst.c464
1 files changed, 464 insertions, 0 deletions
diff --git a/libm/double/mtst.c b/libm/double/mtst.c
new file mode 100644
index 000000000..2559d2340
--- /dev/null
+++ b/libm/double/mtst.c
@@ -0,0 +1,464 @@
+/* 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 <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#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; itst<NTESTS; itst++ )
+{
+d = &defs[itst];
+k = 0;
+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<ntr; 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:
+ 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);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-17 11:35:37 +0000