user/mathtest/mtstf.c

   1 /*   mtst.c
   2  Consistency tests for math functions.
   3  To get strict rounding rules on a 386 or 68000 computer,
   4  define SETPREC to 1.
   5
   6  With NTRIALS=10000, the following are typical results for
   7  IEEE double precision arithmetic.
   8
   9 Consistency test of math functions.
  10 Max and rms relative errors for 10000 random arguments.
  11 x =   cbrt(   cube(x) ):  max = 0.00E+00   rms = 0.00E+00
  12 x =   atan(    tan(x) ):  max = 2.21E-16   rms = 3.27E-17
  13 x =    sin(   asin(x) ):  max = 2.13E-16   rms = 2.95E-17
  14 x =   sqrt( square(x) ):  max = 0.00E+00   rms = 0.00E+00
  15 x =    log(    exp(x) ):  max = 1.11E-16 A rms = 4.35E-18 A
  16 x =   tanh(  atanh(x) ):  max = 2.22E-16   rms = 2.43E-17
  17 x =  asinh(   sinh(x) ):  max = 2.05E-16   rms = 3.49E-18
  18 x =  acosh(   cosh(x) ):  max = 1.43E-15 A rms = 1.54E-17 A
  19 x =  log10(  exp10(x) ):  max = 5.55E-17 A rms = 1.27E-18 A
  20 x = pow( pow(x,a),1/a ):  max = 7.60E-14   rms = 1.05E-15
  21 x =    cos(   acos(x) ):  max = 2.22E-16 A rms = 6.90E-17 A
  22 */
  23
  24 /*
  25 Cephes Math Library Release 2.1:  December, 1988
  26 Copyright 1984, 1987, 1988 by Stephen L. Moshier
  27 */
  28
  29
  30 #include "mconf.h"
  31 #include "math.h"
  32
  33 #define SETPREC 1
  34 #define NTRIALS 10000
  35 #define STRTST 0
  36
  37 #define WTRIALS (NTRIALS/5)
  38
  39 #ifndef ANSIC
  40 float sqrtf(), cbrtf(), expf(), logf();
  41 float exp10f(), log10f(), tanf(), atanf();
  42 float sinf(), asinf(), cosf(), acosf(), powf();
  43 float tanhf(), atanhf(), sinhf(), asinhf(), coshf(), acoshf();
  44 #endif
  45
  46 #define fabsf(x) ((x) < 0 ? -(x) : (x))
  47
  48 #if SETPREC
  49 int sprec();
  50 #endif
  51
  52 int drand();
  53 void exit();
  54 int printf();
  55
  56
  57 /* Provide inverses for square root and cube root: */
  58 #ifdef ANSIC
  59 float square(float x)
  60 #else
  61 float square(x)
  62 float x;
  63 #endif
  64 {
  65 return( x * x );
  66 }
  67
  68 #ifdef ANSIC
  69 float cube(float x)
  70 #else
  71 float cube(x)
  72 float x;
  73 #endif
  74 {
  75 return( x * x * x );
  76 }
  77
  78 /* lookup table for each function */
  79 struct oneargument
  80   {
  81     char *nam1;         /* the function */
  82 #if ANSIC
  83     float (*name) (float);
  84 #else
  85     float (*name) ();
  86 #endif
  87     char *nam2;         /* its inverse  */
  88 #if ANSIC
  89     float (*inv )(float);
  90 #else
  91     float (*inv )();
  92 #endif
  93     int tstyp;          /* type code of the function */
  94     long ctrl;          /* relative error flag */
  95     float arg1w;        /* width of domain for 1st arg */
  96     float arg1l;        /* lower bound domain 1st arg */
  97     long arg1f;         /* flags, e.g. integer arg */
  98   };
  99
 100 struct twoarguments
 101   {
 102     char *nam1;         /* the function */
 103 #if ANSIC
 104     float (*name) (float, float);
 105 #else
 106     float (*name) ();
 107 #endif
 108     char *nam2;         /* its inverse  */
 109 #if ANSIC
 110     float (*inv )(float, float);
 111 #else
 112     float (*inv )();
 113 #endif
 114     int tstyp;          /* type code of the function */
 115     long ctrl;          /* relative error flag */
 116     float arg1w;        /* width of domain for 1st arg */
 117     float arg1l;        /* lower bound domain 1st arg */
 118     long arg1f;         /* flags, e.g. integer arg */
 119     float arg2w;        /* same info for args 2, 3, 4 */
 120     float arg2l;
 121     long arg2f;
 122   };
 123
 124 /* def.ctrl bits: */
 125 #define RELERR 1
 126
 127 /* fundef.tstyp  test types: */
 128 #define POWER 1
 129 #define ELLIP 2
 130 #define GAMMA 3
 131 #define WRONK1 4
 132 #define WRONK2 5
 133 #define WRONK3 6
 134
 135 /* fundef.argNf  argument flag bits: */
 136 #define INT 2
 137 #define EXPSCAL 4
 138
 139 extern float MINLOGF;
 140 extern float MAXLOGF;
 141 extern float PIF;
 142 extern float PIO2F;
 143 /*
 144 define MINLOGF -170.0
 145 define MAXLOGF +170.0
 146 define PIF 3.14159265358979323846
 147 define PIO2F 1.570796326794896619
 148 */
 149
 150 #define N1TESTS 10
 151 struct oneargument defs1arg[N1TESTS] = {
 152 {"  cube",    cube,   "  cbrt",   cbrtf, 0, 1, 2002.0, -1001.0, 0},
 153 {"   tan",    tanf,   "  atan",   atanf, 0, 1,    0.0,     0.0,  0},
 154 {"  asin",   asinf,   "   sin",    sinf, 0, 1,   2.0,      -1.0,  0},
 155 {"square",  square,   "  sqrt",   sqrtf, 0, 1, 87.0,    -43.5, EXPSCAL},
 156 {"   exp",    expf,   "   log",    logf, 0, 0, 174.0,    -87.0,  0},
 157 {" atanh",  atanhf,   "  tanh",   tanhf, 0, 1,    2.0,    -1.0,  0},
 158 {"  sinh",   sinhf,   " asinh",  asinhf, 0, 1, 174.0,   0.0,  0},
 159 {"  cosh",   coshf,   " acosh",  acoshf, 0, 0, 174.0,      0.0,  0},
 160 {" exp10",  exp10f,   " log10",  log10f, 0, 0,  76.0,    -38.0,  0},
 161 {"  acos",   acosf,   "   cos",    cosf, 0, 0,   2.0,      -1.0,  0},
 162 };
 163
 164 #define N2TESTS 1
 165 struct twoarguments defs2arg[N2TESTS] = {
 166 {"pow",       powf,      "pow",    powf, POWER, 1, 20.0, 0.01,   0,
 167 40.0, -20.0, 0},
 168 };
 169
 170 static char *headrs[] = {
 171 "x = %s( %s(x) ): ",
 172 "x = %s( %s(x,a),1/a ): ",      /* power */
 173 "Legendre %s, %s: ",            /* ellip */
 174 "%s(x) = log(%s(x)): ",         /* gamma */
 175 "Wronksian of %s, %s: ",
 176 "Wronksian of %s, %s: ",
 177 "Wronksian of %s, %s: "
 178 };
 179
 180 static float yy1;
 181 static float y2;
 182 static float y3;
 183 static float y4;
 184 static float a;
 185 static float x;
 186 static float y;
 187 static float z;
 188 static float  e;
 189 static float max;
 190 static float rmsa;
 191 static float rms;
 192 static float ave;
 193 static double doublea;
 194
 195 int main()
 196 {
 197 #if ANSIC
 198 float (*fun1 )(float);
 199 float (*ifun1 )(float);
 200 float (*fun2 )(float, float);
 201 float (*ifun2 )(float, float);
 202 #else
 203 float (*fun1 )();
 204 float (*ifun1 )();
 205 float (*fun2 )();
 206 float (*ifun2 )();
 207 #endif
 208 char *nam1, *nam2;
 209 int tstyp, nargs;
 210 long arg1f, arg2f, ctrl;
 211 float arg1l, arg2l, arg1w, arg2w;
 212 int i, k, itst, ntsts, iargs;
 213 int m, ntr;
 214
 215 #if SETPREC
 216 sprec();  /* set coprocessor precision */
 217 #endif
 218 ntr = NTRIALS;
 219 printf( "Consistency test of math functions.\n" );
 220 printf( "Max and rms relative errors for %d random arguments.\n",
 221         ntr );
 222
 223 /* Initialize machine dependent parameters: */
 224
 225 defs1arg[1].arg1w = PIF;
 226 defs1arg[1].arg1l = -PIF/2.0;
 227
 228 /* Microsoft C has trouble with denormal numbers. */
 229 #if 0
 230 defs[3].arg1w = MAXLOGF;
 231 defs[3].arg1l = -MAXLOGF/2.0F;
 232 defs[4].arg1w = 2*MAXLOGF;
 233 defs[4].arg1l = -MAXLOGF;
 234 #endif
 235 defs1arg[6].arg1w = 2.0F*MAXLOGF;
 236 defs1arg[6].arg1l = -MAXLOGF;
 237 defs1arg[7].arg1w = MAXLOGF;
 238 defs1arg[7].arg1l = 0.0;
 239
 240
 241 /* Outer outer loop, on number of function arguments.  */
 242 for( iargs=1; iargs <=2; iargs++)
 243   {
 244     switch (iargs)
 245       {
 246       case 2:
 247         ntsts = N2TESTS;
 248         break;
 249       default:
 250         ntsts = N1TESTS;
 251       }
 252 /* Outer loop, on the test number: */
 253 for( itst=STRTST; itst<ntsts; itst++ )
 254 {
 255   switch (iargs)
 256     {
 257     case 2:
 258       tstyp = defs2arg[itst].tstyp;
 259       fun2 = defs2arg[itst].name;
 260       ifun2 = defs2arg[itst].inv;
 261       nam1 = defs2arg[itst].nam1;
 262       nam2 = defs2arg[itst].nam2;
 263       arg1w = defs2arg[itst].arg1w;
 264       arg1l = defs2arg[itst].arg1l;
 265       arg1f = defs2arg[itst].arg1f;
 266       arg2w = defs2arg[itst].arg2w;
 267       arg2l = defs2arg[itst].arg2l;
 268       arg2f = defs2arg[itst].arg2f;
 269       ctrl = defs2arg[itst].ctrl;
 270       nargs = 2;
 271       break;
 272     default:
 273       tstyp = defs1arg[itst].tstyp;
 274       fun1 = defs1arg[itst].name;
 275       ifun1 = defs1arg[itst].inv;
 276       nam1 = defs1arg[itst].nam1;
 277       nam2 = defs1arg[itst].nam2;
 278       arg1w = defs1arg[itst].arg1w;
 279       arg1l = defs1arg[itst].arg1l;
 280       arg1f = defs1arg[itst].arg1f;
 281       ctrl = defs1arg[itst].ctrl;
 282       nargs = 1;
 283     }
 284 k = 0;
 285 m = 0;
 286 max = 0.0F;
 287 rmsa = 0.0F;
 288 ave = 0.0F;
 289
 290 /* Absolute error criterion starts with gamma function
 291  * (put all such at end of table)
 292  */
 293 if( tstyp == GAMMA )
 294         printf( "Absolute error criterion (but relative if >1):\n" );
 295
 296 /* Smaller number of trials for Wronksians
 297  * (put them at end of list)
 298  */
 299 if( tstyp == WRONK1 )
 300         {
 301         ntr = WTRIALS;
 302         printf( "Absolute error and only %d trials:\n", ntr );
 303         }
 304
 305 printf( headrs[tstyp], nam2, nam1 );
 306
 307 for( i=0; i<ntr; i++ )
 308 {
 309 m++;
 310
 311 /* make random number(s) in desired range(s) */
 312 switch( nargs )
 313 {
 314
 315 default:
 316 goto illegn;
 317
 318 case 2:
 319 drand( &doublea );
 320 a = arg2w *  ( doublea - 1.0 )  +  arg2l;
 321 if( arg2f & EXPSCAL )
 322         {
 323         a = expf(a);
 324         drand( &doublea );
 325         y2 = doublea;
 326         a -= 1.0e-13 * a * y2;
 327         }
 328 if( arg2f & INT )
 329         {
 330         k = a + 0.25;
 331         a = k;
 332         }
 333
 334 case 1:
 335 drand( &doublea );
 336 x = arg1w *  ( doublea - 1.0 )  +  arg1l;
 337 if( arg1f & EXPSCAL )
 338         {
 339         x = expf(x);
 340         drand( &doublea );
 341         a = doublea;
 342         x += 1.0e-13F * x * a;
 343         }
 344 }
 345
 346
 347 /* compute function under test */
 348 switch( nargs )
 349         {
 350         case 2:
 351         if( arg2f & INT )
 352                 {
 353                 switch( tstyp )
 354                         {
 355                         case WRONK1:
 356                         yy1 = (*fun2)( k, x ); /* jn */
 357                         y2 = (*fun2)( k+1, x );
 358                         y3 = (*ifun2)( k, x ); /* yn */
 359                         y4 = (*ifun2)( k+1, x );
 360                         break;
 361
 362                         case WRONK2:
 363                         yy1 = (*fun2)( a, x ); /* iv */
 364                         y2 = (*fun2)( a+1.0F, x );
 365                         y3 = (*ifun2)( k, x ); /* kn */
 366                         y4 = (*ifun2)( k+1, x );
 367                         break;
 368
 369                         default:
 370                         z = (*fun2)( k, x );
 371                         y = (*ifun2)( k, z );
 372                         }
 373                 }
 374         else
 375                 {
 376                 if( tstyp == POWER )
 377                         {
 378                         z = (*fun2)( x, a );
 379                         y = (*ifun2)( z, 1.0F/a );
 380                         }
 381                 else
 382                         {
 383                         z = (*fun2)( a, x );
 384                         y = (*ifun2)( a, z );
 385                         }
 386                 }
 387         break;
 388
 389         case 1:
 390         switch( tstyp )
 391                 {
 392                 case ELLIP:
 393                 yy1 = ( *(fun1) )(x);
 394                 y2 = ( *(fun1) )(1.0F-x);
 395                 y3 = ( *(ifun1) )(x);
 396                 y4 = ( *(ifun1) )(1.0F-x);
 397                 break;
 398
 399 #if 0
 400                 case GAMMA:
 401                 y = lgam(x);
 402                 x = log( gamma(x) );
 403                 break;
 404 #endif
 405                 default:
 406                 z = ( *(fun1) )(x);
 407                 y = ( *(ifun1) )(z);
 408                 }
 409         break;
 410
 411
 412         default:
 413 illegn:
 414         printf( "Illegal nargs= %d", nargs );
 415         exit(1);
 416         }
 417
 418 switch( tstyp )
 419         {
 420         case WRONK1:
 421         e = (y2*y3 - yy1*y4) - 2.0F/(PIF*x); /* Jn, Yn */
 422         break;
 423
 424         case WRONK2:
 425         e = (y2*y3 + yy1*y4) - 1.0F/x; /* In, Kn */
 426         break;
 427
 428         case ELLIP:
 429         e = (yy1-y3)*y4 + y3*y2 - PIO2F;
 430         break;
 431
 432         default:
 433         e = y - x;
 434         break;
 435         }
 436
 437 if( ctrl & RELERR )
 438         e /= x;
 439 else
 440         {
 441         if( fabsf(x) > 1.0F )
 442                 e /= x;
 443         }
 444
 445 ave += e;
 446 /* absolute value of error */
 447 if( e < 0 )
 448         e = -e;
 449
 450 /* peak detect the error */
 451 if( e > max )
 452         {
 453         max = e;
 454
 455         if( e > 1.0e-3F )
 456                 {
 457                 printf("x %.6E z %.6E y %.6E max %.4E\n",
 458                  x, z, y, max);
 459                 if( tstyp == POWER )
 460                         {
 461                         printf( "a %.6E\n", a );
 462                         }
 463                 if( tstyp >= WRONK1 )
 464                         {
 465                 printf( "yy1 %.4E y2 %.4E y3 %.4E y4 %.4E k %d x %.4E\n",
 466                  yy1, y2, y3, y4, k, x );
 467                         }
 468                 }
 469
 470 /*
 471         printf("%.8E %.8E %.4E %6ld \n", x, y, max, n);
 472         printf("%d %.8E %.8E %.4E %6ld \n", k, x, y, max, n);
 473         printf("%.6E %.6E %.6E %.4E %6ld \n", a, x, y, max, n);
 474         printf("%.6E %.6E %.6E %.6E %.4E %6ld \n", a, b, x, y, max, n);
 475         printf("%.4E %.4E %.4E %.4E %.4E %.4E %6ld \n",
 476                 a, b, c, x, y, max, n);
 477 */
 478         }
 479
 480 /* accumulate rms error */
 481 e *= 1.0e7F;    /* adjust range */
 482 rmsa += e * e;  /* accumulate the square of the error */
 483 }
 484
 485 /* report after NTRIALS trials */
 486 rms = 1.0e-7F * sqrtf( rmsa/m );
 487 if(ctrl & RELERR)
 488         printf(" max = %.10E   rms = %.10E\n", max, rms );
 489 else
 490         printf(" max = %.10E A rms = %.10E A\n", max, rms );
 491 } /* loop on itst */
 492 } /* loop on number of args */
 493 return 0;
 494 }