test/math/ieetst.c

   1 /* Floating point to ASCII input and output string test program.\r
   2  *\r
   3  * Numbers in the native machine data structure are converted\r
   4  * to e type, then to and from decimal ASCII strings.  Native\r
   5  * printf() and scanf() functions are also used to produce\r
   6  * and read strings.  The resulting e type binary values\r
   7  * are compared, with diagnostic printouts of any discrepancies.\r
   8  *\r
   9  * Steve Moshier, 16 Dec 88\r
  10  * last revision: 16 May 92\r
  11  */\r
  12 \r
  13 #include "ehead.h"\r
  14 #include "mconf.h"\r
  15 \r
  16 /* Include tests of 80-bit long double precision: */\r
  17 #define LDOUBLE 0\r
  18 /* Abort subtest after getting this many errors: */\r
  19 #define MAXERR 5\r
  20 /* Number of random arguments to try (set as large as you have\r
  21  * patience for): */\r
  22 #define NRAND 100\r
  23 /* Perform internal consistency test: */\r
  24 #define CHKINTERNAL 0\r
  25 \r
  26 static unsigned short fullp[NE], rounded[NE];\r
  27 float prec24, sprec24, ssprec24;\r
  28 double prec53, sprec53, ssprec53;\r
  29 #if LDOUBLE\r
  30 long double prec64, sprec64, ssprec64;\r
  31 #endif\r
  32 \r
  33 static unsigned short rprint[NE], rscan[NE];\r
  34 static unsigned short q1[NE], q2[NE], q5[NE];\r
  35 static unsigned short e1[NE], e2[NE], e3[NE];\r
  36 static double d1, d2;\r
  37 static int errprint = 0;\r
  38 static int errscan = 0;\r
  39 static int identerr = 0;\r
  40 static int errtot = 0;\r
  41 static int count = 0;\r
  42 static char str0[80], str1[80], str2[80], str3[80];\r
  43 static unsigned short eten[NE], maxm[NE];\r
  44 \r
  45 int m, n, k2, mprec, SPREC;\r
  46 \r
  47 char *Ten = "10.0";\r
  48 char tformat[10];\r
  49 char *format24 = "%.8e";\r
  50 #ifdef DEC\r
  51 char *format53 = "%.17e";\r
  52 #else\r
  53 char *format53 = "%.16e";\r
  54 #endif\r
  55 char *fformat24 = "%e";\r
  56 char *fformat53 = "%le";\r
  57 char *pct = "%";\r
  58 char *quo = "\042";\r
  59 #if LDOUBLE\r
  60 char *format64 = "%.20Le";\r
  61 char *fformat64 = "%Le";\r
  62 #endif\r
  63 char *format;\r
  64 char *fformat;\r
  65 char *toomany = "Too many errors; aborting this test.\n";\r
  66 \r
  67 static int mnrflag;\r
  68 static int etrflag;\r
  69 void chkit(), printerr(), mnrand(), etrand(), shownoncrit();\r
  70 void chkid(), pvec();\r
  71 \r
  72 main()\r
  73 {\r
  74 int i, iprec;\r
  75 \r
  76 printf( "Steve Moshier's printf/scanf tester, version 0.2.\n\n" );\r
  77 #ifdef DEC\r
  78  /* DEC PDP-11/VAX single precision not yet implemented */\r
  79 for( iprec = 1; iprec<2; iprec++ )\r
  80 #else\r
  81 for( iprec = 0; iprec<3; iprec++ )\r
  82 #endif\r
  83         {\r
  84         errscan = 0;\r
  85         identerr = 0;\r
  86         errprint = 0;\r
  87         eclear( rprint );\r
  88         eclear( rscan );\r
  89 \r
  90 switch( iprec )\r
  91         {\r
  92         case 0:\r
  93                 SPREC = 8; /* # digits after the decimal point */\r
  94                 mprec = 24; /* # bits in the significand */\r
  95                 m = 9; /* max # decimal digits for correct rounding */\r
  96                 n = 13; /* max power of ten for correct rounding */\r
  97                 k2 = -125; /* underflow beyond 2^-k2 */\r
  98                 format = format24; /* printf format string */\r
  99                 fformat = fformat24; /* scanf format string */\r
 100                 mnrflag = 1; /* sets interval for random numbers */\r
 101                 etrflag = 1;\r
 102                 printf( "Testing FLOAT precision.\n" );\r
 103                 break;\r
 104 \r
 105         case 1:\r
 106 #ifdef DEC\r
 107                 SPREC = 17;\r
 108                 mprec = 56;\r
 109                 m = 17;\r
 110                 n = 27;\r
 111                 k2 = -125;\r
 112                 format = format53;\r
 113                 fformat = fformat53;\r
 114                 mnrflag = 2;\r
 115                 etrflag = 1;\r
 116                 printf( "Testing DEC DOUBLE precision.\n" );\r
 117                 break;\r
 118 #else\r
 119                 SPREC = 16;\r
 120                 mprec = 53;\r
 121                 m = 17;\r
 122                 n = 27;\r
 123                 k2 = -1021;\r
 124                 format = format53;\r
 125                 fformat = fformat53;\r
 126                 mnrflag = 2;\r
 127                 etrflag = 2;\r
 128                 printf( "Testing DOUBLE precision.\n" );\r
 129                 break;\r
 130 #endif\r
 131         case 2:\r
 132 #if LDOUBLE\r
 133                 SPREC = 20;\r
 134                 mprec = 64;\r
 135                 m = 20;\r
 136                 n = 34;\r
 137                 k2 = -16382;\r
 138                 format = format64;\r
 139                 fformat = fformat64;\r
 140                 mnrflag = 3;\r
 141                 etrflag = 3;\r
 142                 printf( "Testing LONG DOUBLE precision.\n" );\r
 143                 break;\r
 144 #else\r
 145                 goto nodenorm;\r
 146 #endif\r
 147         }\r
 148 \r
 149         asctoe( Ten, eten );\r
 150 /* 10^m - 1 */\r
 151         d2 = m;\r
 152         e53toe( &d2, e1 );\r
 153         epow( eten, e1, maxm );\r
 154         esub( eone, maxm, maxm );\r
 155 \r
 156 /* test 1 */\r
 157         printf( "1. Checking 10^n - 1 for n = %d to %d.\n", -m, m );\r
 158         emov( eone, q5 );\r
 159         for( count=0; count<=m; count++ )\r
 160                 {\r
 161                 esub( eone, q5, fullp );\r
 162                 chkit( 1 );\r
 163                 ediv( q5, eone, q2 );\r
 164                 esub( eone, q2, fullp );\r
 165                 chkit( 1 );\r
 166                 emul( eten, q5, q5 );\r
 167                 if( errtot >= MAXERR )\r
 168                         {\r
 169                         printf( "%s", toomany );\r
 170                         goto end1;\r
 171                         }\r
 172                 }\r
 173 end1:\r
 174         printerr();\r
 175 \r
 176 \r
 177 /* test 2 */\r
 178         printf( "2. Checking powers of 10 from 10^-%d to 10^%d.\n", n, n );\r
 179         emov( eone, q5 );\r
 180         for( count=0; count<=n; count++ )\r
 181                 {\r
 182                 emov( q5, fullp );\r
 183                 chkit( 2 );\r
 184                 ediv( q5, eone, fullp );\r
 185                 chkit( 2 );\r
 186                 emul( eten, q5, q5 );\r
 187                 if( errtot >= MAXERR )\r
 188                         {\r
 189                         printf( "%s", toomany );\r
 190                         goto end2;\r
 191                         }\r
 192                 }\r
 193 end2:\r
 194         printerr();\r
 195 \r
 196 /* test 3 */\r
 197         printf( "3. Checking (10^%d-1)*10^n from n = -%d to %d.\n", m, n, n );\r
 198         emov( eone, q5 );\r
 199         for( count= -n; count<=n; count++ )\r
 200                 {\r
 201                 emul( maxm, q5, fullp );\r
 202                 chkit( 3 );\r
 203                 emov( q5, fullp );\r
 204                 ediv( fullp, eone, fullp );\r
 205                 emul( maxm, fullp, fullp );\r
 206                 chkit( 3 );\r
 207                 emul( eten, q5, q5 );\r
 208                 if( errtot >= MAXERR )\r
 209                         {\r
 210                         printf( "%s", toomany );\r
 211                         goto end3;\r
 212                         }\r
 213                 }\r
 214 end3:\r
 215         printerr();\r
 216 \r
 217 \r
 218 \r
 219 /* test 4 */\r
 220         printf( "4. Checking powers of 2 from 2^-24 to 2^+56.\n" );\r
 221         d1 = -24.0;\r
 222         e53toe( &d1, q1 );\r
 223         epow( etwo, q1, q5 );\r
 224 \r
 225         for( count = -24; count <= 56; count++ )\r
 226                 {\r
 227                 emov( q5, fullp );\r
 228                 chkit( 4 );\r
 229                 emul( etwo, q5, q5 );\r
 230                 if( errtot >= MAXERR )\r
 231                         {\r
 232                         printf( "%s", toomany );\r
 233                         goto end4;\r
 234                         }\r
 235                 }\r
 236 end4:\r
 237         printerr();\r
 238 \r
 239 \r
 240 /* test 5 */\r
 241         printf( "5. Checking 2^n - 1 for n = 0 to %d.\n", mprec );\r
 242         emov( eone, q5 );\r
 243         for( count=0; count<=mprec; count++ )\r
 244                 {\r
 245                 esub( eone, q5, fullp );\r
 246                 chkit( 5 );\r
 247                 emul( etwo, q5, q5 );\r
 248                 if( errtot >= MAXERR )\r
 249                         {\r
 250                         printf( "%s", toomany );\r
 251                         goto end5;\r
 252                         }\r
 253                 }\r
 254 end5:\r
 255         printerr();\r
 256 \r
 257 /* test 6 */\r
 258         printf( "6. Checking 2^n + 1 for n = 0 to %d.\n", mprec );\r
 259         emov( eone, q5 );\r
 260         for( count=0; count<=mprec; count++ )\r
 261                 {\r
 262                 eadd( eone, q5, fullp );\r
 263                 chkit( 6 );\r
 264                 emul( etwo, q5, q5 );\r
 265                 if( errtot >= MAXERR )\r
 266                         {\r
 267                         printf( "%s", toomany );\r
 268                         goto end6;\r
 269                         }\r
 270                 }\r
 271 end6:\r
 272         printerr();\r
 273 \r
 274 /* test 7 */\r
 275         printf(\r
 276          "7. Checking %d values M * 10^N with random integer M and N,\n",\r
 277          NRAND );\r
 278         printf("  1 <= M <= 10^%d - 1  and  -%d <= N <= +%d.\n", m, n, n );\r
 279         for( i=0; i<NRAND; i++ )\r
 280                 {\r
 281                 mnrand( fullp );\r
 282                 chkit( 7 );\r
 283                 if( errtot >= MAXERR )\r
 284                         {\r
 285                         printf( "%s", toomany );\r
 286                         goto end7;\r
 287                         }\r
 288                 }\r
 289 end7:\r
 290         printerr();\r
 291 \r
 292 /* test 8 */\r
 293         printf("8. Checking critical rounding cases.\n" );\r
 294         for( i=0; i<20; i++ )\r
 295                 {\r
 296                 mnrand( fullp );\r
 297                 eabs( fullp );\r
 298                 if( ecmp( fullp, eone ) < 0 )\r
 299                         ediv( fullp, eone, fullp );\r
 300                 efloor( fullp, fullp );\r
 301                 eadd( ehalf, fullp, fullp );\r
 302                 chkit( 8 );\r
 303                 if( errtot >= MAXERR )\r
 304                         {\r
 305                         printf( "%s", toomany );\r
 306                         goto end8;\r
 307                         }\r
 308                 }\r
 309 end8:\r
 310         printerr();\r
 311 \r
 312 \r
 313 \r
 314 /* test 9 */\r
 315         printf("9. Testing on %d random non-denormal values.\n", NRAND );\r
 316         for( i=0; i<NRAND; i++ )\r
 317                 {\r
 318                 etrand( fullp );\r
 319                 chkit( 9 );\r
 320                 }\r
 321         printerr();\r
 322         shownoncrit();\r
 323 \r
 324 /* test 10 */\r
 325         printf(\r
 326         "Do you want to check denormal numbers in this precision ? (y/n) " );\r
 327         gets( str0 );\r
 328         if( str0[0] != 'y' )\r
 329                 goto nodenorm;\r
 330 \r
 331         printf( "10. Checking denormal numbers.\n" );\r
 332 \r
 333 /* Form 2^-starting power */\r
 334         d1 = k2;\r
 335         e53toe( &d1, q1 );\r
 336         epow( etwo, q1, e1 );\r
 337 \r
 338 /* Find 2^-mprec less than starting power */\r
 339         d1 = -mprec + 4;\r
 340         e53toe( &d1, q1 );\r
 341         epow( etwo, q1, e3 );\r
 342         emul( e1, e3, e3 );\r
 343         emov( e3, e2 );\r
 344         ediv( etwo, e2, e2 );\r
 345 \r
 346         while( ecmp(e1,e2) != 0 )\r
 347                 {\r
 348                 eadd( e1, e2, fullp );\r
 349                 switch( mprec )\r
 350                         {\r
 351 #if LDOUBLE\r
 352                         case 64:\r
 353                         etoe64( e1, &sprec64 );\r
 354                         e64toe( &sprec64, q1 );\r
 355                         etoe64( fullp, &prec64 );\r
 356                         e64toe( &prec64, q2 );\r
 357                         break;\r
 358 #endif\r
 359 #ifdef DEC\r
 360                         case 56:\r
 361 #endif\r
 362                         case 53:\r
 363                         etoe53( e1, &sprec53 );\r
 364                         e53toe( &sprec53, q1 );\r
 365                         etoe53( fullp, &prec53 );\r
 366                         e53toe( &prec53, q2 );\r
 367                         break;\r
 368 \r
 369                         case 24:\r
 370                         etoe24( e1, &sprec24 );\r
 371                         e24toe( &sprec24, q1 );\r
 372                         etoe24( fullp, &prec24 );\r
 373                         e24toe( &prec24, q2 );\r
 374                         break;\r
 375                         }\r
 376                 if( ecmp( q2, ezero ) == 0 )\r
 377                         goto maxden;\r
 378                 chkit(10);\r
 379                 if( ecmp(q1,q2) == 0 )\r
 380                         {\r
 381                         ediv( etwo, e1, e1 );\r
 382                         emov( e3, e2 );\r
 383                         }\r
 384                 if( errtot >= MAXERR )\r
 385                         {\r
 386                         printf( "%s", toomany );\r
 387                         goto maxden;\r
 388                         }\r
 389                 ediv( etwo, e2, e2 );\r
 390                 }\r
 391 maxden:\r
 392         printerr();\r
 393 nodenorm:\r
 394         printf( "\n" );\r
 395         } /* loop on precision */\r
 396 printf( "End of test.\n" );\r
 397 }\r
 398 \r
 399 #if CHKINTERNAL\r
 400 long double xprec64;\r
 401 double xprec53;\r
 402 float xprec24;\r
 403 \r
 404 /* Check binary -> printf -> scanf -> binary identity\r
 405  * of internal routines\r
 406  */\r
 407 void chkinternal( ref, tst, string )\r
 408 unsigned short ref[], tst[];\r
 409 char *string;\r
 410 {\r
 411 \r
 412 if( ecmp(ref,tst) != 0 )\r
 413         {\r
 414         printf( "internal identity compare error!\n" );\r
 415         chkid( ref, tst, string );\r
 416         }\r
 417 }\r
 418 #endif\r
 419 \r
 420 \r
 421 /* Check binary -> printf -> scanf -> binary identity\r
 422  */\r
 423 void chkid( print, scan, string )\r
 424 unsigned short print[], scan[];\r
 425 char *string;\r
 426 {\r
 427 /* Test printf-scanf identity */\r
 428 if( ecmp( print, scan ) != 0 )\r
 429         {\r
 430         pvec( print, NE );\r
 431         printf( " ->printf-> %s ->scanf->\n", string );\r
 432         pvec( scan, NE );\r
 433         printf( " is not an identity.\n" );\r
 434         ++identerr;\r
 435         }\r
 436 }\r
 437 \r
 438 \r
 439 /* Check scanf result\r
 440  */\r
 441 void chkscan( ref, tst, string )\r
 442 unsigned short ref[], tst[];\r
 443 char *string;\r
 444 {\r
 445 /* Test scanf()  */\r
 446 if( ecmp( ref, tst ) != 0 )\r
 447         {\r
 448         printf( "scanf(%s) -> ", string );\r
 449         pvec( tst, NE );\r
 450         printf( "\n should be    " );\r
 451         pvec( ref, NE );\r
 452         printf( ".\n" );\r
 453         ++errscan;\r
 454         ++errtot;\r
 455         }\r
 456 }\r
 457 \r
 458 \r
 459 /* Test printf() result\r
 460  */\r
 461 void chkprint( ref, tst, string ) \r
 462 unsigned short ref[], tst[];\r
 463 char *string;\r
 464 {\r
 465 if( ecmp(ref, tst) != 0 )\r
 466         {\r
 467         printf( "printf( ");\r
 468         pvec( ref, NE );\r
 469         printf( ") -> %s\n", string );\r
 470         printf( "      = " );\r
 471         pvec( tst, NE );\r
 472         printf( ".\n" );\r
 473         ++errprint;\r
 474         ++errtot;\r
 475         }\r
 476 }\r
 477 \r
 478 \r
 479 /* Print array of n 16-bit shorts\r
 480  */\r
 481 void pvec( x, n )\r
 482 unsigned short x[];\r
 483 int n;\r
 484 {\r
 485 int i;\r
 486 \r
 487 for( i=0; i<n; i++ )\r
 488         {\r
 489         printf( "%04x ", x[i] );\r
 490         }\r
 491 }\r
 492 \r
 493 /* Measure worst case printf rounding error\r
 494  */\r
 495 void cmpprint( ref, tst )\r
 496 unsigned short ref[], tst[];\r
 497 {\r
 498 unsigned short e[NE];\r
 499 \r
 500 if( ecmp( ref, ezero ) != 0 )\r
 501         {\r
 502         esub( ref, tst, e );\r
 503         ediv( ref, e, e );\r
 504         eabs( e );\r
 505         if( ecmp( e, rprint ) > 0 )\r
 506                 emov( e, rprint );\r
 507         }\r
 508 }\r
 509 \r
 510 /* Measure worst case scanf rounding error\r
 511  */\r
 512 void cmpscan( ref, tst )\r
 513 unsigned short ref[], tst[];\r
 514 {\r
 515 unsigned short er[NE];\r
 516 \r
 517 if( ecmp( ref, ezero ) != 0 )\r
 518         {\r
 519         esub( ref, tst, er );\r
 520         ediv( ref, er, er );\r
 521         eabs( er );\r
 522         if( ecmp( er, rscan ) > 0 )\r
 523                 emov( er, rscan );\r
 524         if( ecmp( er, ehalf ) > 0 )\r
 525                 {\r
 526                 etoasc( tst, str1, 21 );\r
 527                 printf( "Bad error: scanf(%s) = %s !\n", str0, str1 );\r
 528                 }\r
 529         }\r
 530 }\r
 531 \r
 532 /* Check rounded-down decimal string output of printf\r
 533  */\r
 534 void cmptrunc( ref, tst )\r
 535 unsigned short ref[], tst[];\r
 536 {\r
 537 if( ecmp( ref, tst ) != 0 )\r
 538         {\r
 539         printf( "printf(%s%s%s, %s) -> %s\n", quo, tformat, quo, str1, str2 );\r
 540         printf( "should be      %s .\n", str3 );\r
 541         errprint += 1;\r
 542         }\r
 543 }\r
 544 \r
 545 \r
 546 void shownoncrit()\r
 547 {\r
 548 \r
 549 etoasc( rprint, str0, 3 );\r
 550 printf( "Maximum relative printf error found = %s .\n", str0 );\r
 551 etoasc( rscan, str0, 3 );\r
 552 printf( "Maximum relative scanf error found = %s .\n", str0 );\r
 553 }\r
 554 \r
 555 \r
 556 \r
 557 /* Produce arguments and call comparison subroutines.\r
 558  */\r
 559 void chkit( testno )\r
 560 int testno;\r
 561 {\r
 562 unsigned short t[NE], u[NE], v[NE];\r
 563 int j;\r
 564 \r
 565 switch( mprec )\r
 566         {\r
 567 #if LDOUBLE\r
 568         case 64:\r
 569                 etoe64( fullp, &prec64 );\r
 570                 e64toe( &prec64, rounded );\r
 571 #if CHKINTERNAL\r
 572                 e64toasc( &prec64, str1, SPREC );\r
 573                 asctoe64( str1, &xprec64 );\r
 574                 e64toe( &xprec64, t );\r
 575                 chkinternal( rounded, t, str1 );\r
 576 #endif\r
 577 /* check printf and scanf */\r
 578                 sprintf( str2, format, prec64 );\r
 579                 sscanf( str2, fformat, &sprec64 );\r
 580                 e64toe( &sprec64, u );\r
 581                 chkid( rounded, u, str2 );\r
 582                 asctoe64( str2, &ssprec64 );\r
 583                 e64toe( &ssprec64, v );\r
 584                 chkscan( v, u, str2 );\r
 585                 chkprint( rounded, v, str2 );\r
 586                 if( testno < 8 )\r
 587                         break;\r
 588 /* rounding error measurement */\r
 589                 etoasc( fullp, str0, 24 );\r
 590                 etoe64( fullp, &ssprec64 );\r
 591                 e64toe( &ssprec64, u );\r
 592                 sprintf( str2, format, ssprec64 );\r
 593                 asctoe( str2, t );\r
 594                 cmpprint( u, t );\r
 595                 sscanf( str0, fformat, &sprec64 );\r
 596                 e64toe( &sprec64, t );\r
 597                 cmpscan( fullp, t );\r
 598                 if( testno < 8 )\r
 599                         break;\r
 600 /* strings rounded to less than maximum precision */\r
 601                 e64toasc( &ssprec64, str1, 24 );\r
 602                 for( j=SPREC-1; j>0; j-- )              \r
 603                         {\r
 604                         e64toasc( &ssprec64, str3, j );\r
 605                         asctoe( str3, v );\r
 606                         sprintf( tformat, "%s.%dLe", pct, j );\r
 607                         sprintf( str2, tformat, ssprec64 );\r
 608                         asctoe( str2, t );\r
 609                         cmptrunc( v, t );\r
 610                         }\r
 611                 break;\r
 612 #endif\r
 613 #ifdef DEC\r
 614         case 56:\r
 615 #endif\r
 616         case 53:\r
 617                 etoe53( fullp, &prec53 );\r
 618                 e53toe( &prec53, rounded );\r
 619 #if CHKINTERNAL\r
 620                 e53toasc( &prec53, str1, SPREC );\r
 621                 asctoe53( str1, &xprec53 );\r
 622                 e53toe( &xprec53, t );\r
 623                 chkinternal( rounded, t, str1 );\r
 624 #endif\r
 625                 sprintf( str2, format, prec53 );\r
 626                 sscanf( str2, fformat, &sprec53 );\r
 627                 e53toe( &sprec53, u );\r
 628                 chkid( rounded, u, str2 );\r
 629                 asctoe53( str2, &ssprec53 );\r
 630                 e53toe( &ssprec53, v );\r
 631                 chkscan( v, u, str2 );\r
 632                 chkprint( rounded, v, str2 );\r
 633                 if( testno < 8 )\r
 634                         break;\r
 635 /* rounding error measurement */\r
 636                 etoasc( fullp, str0, 24 );\r
 637                 etoe53( fullp, &ssprec53 );\r
 638                 e53toe( &ssprec53, u );\r
 639                 sprintf( str2, format, ssprec53 );\r
 640                 asctoe( str2, t );\r
 641                 cmpprint( u, t );\r
 642                 sscanf( str0, fformat, &sprec53 );\r
 643                 e53toe( &sprec53, t );\r
 644                 cmpscan( fullp, t );\r
 645                 if( testno < 8 )\r
 646                         break;\r
 647                 e53toasc( &ssprec53, str1, 24 );\r
 648                 for( j=SPREC-1; j>0; j-- )              \r
 649                         {\r
 650                         e53toasc( &ssprec53, str3, j );\r
 651                         asctoe( str3, v );\r
 652                         sprintf( tformat, "%s.%de", pct, j );\r
 653                         sprintf( str2, tformat, ssprec53 );\r
 654                         asctoe( str2, t );\r
 655                         cmptrunc( v, t );\r
 656                         }\r
 657                 break;\r
 658 \r
 659         case 24:\r
 660                 etoe24( fullp, &prec24 );\r
 661                 e24toe( &prec24, rounded );\r
 662 #if CHKINTERNAL\r
 663                 e24toasc( &prec24, str1, SPREC );\r
 664                 asctoe24( str1, &xprec24 );\r
 665                 e24toe( &xprec24, t );\r
 666                 chkinternal( rounded, t, str1 );\r
 667 #endif\r
 668                 sprintf( str2, format, prec24 );\r
 669                 sscanf( str2, fformat, &sprec24 );\r
 670                 e24toe( &sprec24, u );\r
 671                 chkid( rounded, u, str2 );\r
 672                 asctoe24( str2, &ssprec24 );\r
 673                 e24toe( &ssprec24, v );\r
 674                 chkscan( v, u, str2 );\r
 675                 chkprint( rounded, v, str2 );\r
 676                 if( testno < 8 )\r
 677                         break;\r
 678 /* rounding error measurement */\r
 679                 etoasc( fullp, str0, 24 );\r
 680                 etoe24( fullp, &ssprec24 );\r
 681                 e24toe( &ssprec24, u );\r
 682                 sprintf( str2, format, ssprec24 );\r
 683                 asctoe( str2, t );\r
 684                 cmpprint( u, t );\r
 685                 sscanf( str0, fformat, &sprec24 );\r
 686                 e24toe( &sprec24, t );\r
 687                 cmpscan( fullp, t );\r
 688 /*\r
 689                 if( testno < 8 )\r
 690                         break;\r
 691 */\r
 692                 e24toasc( &ssprec24, str1, 24 );\r
 693                 for( j=SPREC-1; j>0; j-- )              \r
 694                         {\r
 695                         e24toasc( &ssprec24, str3, j );\r
 696                         asctoe( str3, v );\r
 697                         sprintf( tformat, "%s.%de", pct, j );\r
 698                         sprintf( str2, tformat, ssprec24 );\r
 699                         asctoe( str2, t );\r
 700                         cmptrunc( v, t );\r
 701                         }\r
 702                 break;\r
 703         }\r
 704 }\r
 705 \r
 706 \r
 707 void printerr()\r
 708 {\r
 709 if( (errscan == 0) && (identerr == 0) && (errprint == 0) )\r
 710         printf( "No errors found.\n" );\r
 711 else\r
 712         {\r
 713         printf( "%d binary -> decimal errors found.\n", errprint );\r
 714         printf( "%d decimal -> binary errors found.\n", errscan );\r
 715         }\r
 716 errscan = 0;    /* reset for next test */\r
 717 identerr = 0;\r
 718 errprint = 0;\r
 719 errtot = 0;\r
 720 }\r
 721 \r
 722 \r
 723 /* Random number generator\r
 724  * in the range M * 10^N, where 1 <= M <= 10^17 - 1\r
 725  * and -27 <= N <= +27.  Test values of M are logarithmically distributed\r
 726  * random integers; test values of N are uniformly distributed random integers.\r
 727  */\r
 728 \r
 729 static char *fwidth = "1.036163291797320557783096e1"; /* log(sqrt(10^9-1)) */\r
 730 static char *dwidth = "1.957197329044938830915E1"; /* log(sqrt(10^17-1)) */\r
 731 static char *ldwidth = "2.302585092994045684017491e1"; /* log(sqrt(10^20-1)) */\r
 732 \r
 733 static char *a13 = "13.0";\r
 734 static char *a27 = "27.0";\r
 735 static char *a34 = "34.0";\r
 736 static char *a10m13 = "1.0e-13";\r
 737 static unsigned short LOW[ NE ], WIDTH[NE], e27[NE], e10m13[NE];\r
 738 \r
 739 \r
 740 void mnrand( erand )\r
 741 unsigned short erand[];\r
 742 {\r
 743 unsigned short ea[NE], em[NE], en[NE], ex[NE];\r
 744 double x, a;\r
 745 \r
 746 if( mnrflag )\r
 747         {\r
 748         if( mnrflag == 3 )\r
 749                 {\r
 750                 asctoe( ldwidth, WIDTH );\r
 751                 asctoe( a34, e27 );\r
 752                 }\r
 753         if( mnrflag == 2 )\r
 754                 {\r
 755                 asctoe( dwidth, WIDTH );\r
 756                 asctoe( a27, e27 );\r
 757                 }\r
 758         if( mnrflag == 1 )\r
 759                 {\r
 760                 asctoe( fwidth, WIDTH );\r
 761                 asctoe( a13, e27 );\r
 762                 }\r
 763         asctoe( a10m13, e10m13 );\r
 764         mnrflag = 0;\r
 765         }\r
 766 drand( &x );\r
 767 e53toe( &x, ex ); /* x = WIDTH *  ( x - 1.0 )  +  LOW; */\r
 768 esub( eone, ex, ex );\r
 769 emul( WIDTH, ex, ex );\r
 770 eexp( ex, ex );   /* x = exp(x); */\r
 771 \r
 772 drand( &a );\r
 773 e53toe( &a, ea );\r
 774 emul( ea, ex, ea );  /* a = 1.0e-13 * x * a; */\r
 775 emul( e10m13, ea, ea );\r
 776 eabs( ea );\r
 777 eadd( ea, ex, ex );     /* add fuzz */\r
 778 emul( ex, ex, ex );     /* square it, to get range to 10^17 - 1 */\r
 779 efloor( ex, em ); /* this is M */\r
 780 \r
 781 /* Random power of 10 */\r
 782 drand( &a );\r
 783 e53toe( &a, ex );\r
 784 esub( eone, ex, ex ); /* y3 = 54.0 *  ( y3 - 1.0 ) + 0.5; */\r
 785 emul( e27, ex, ex );\r
 786 eadd( ex, ex, ex );\r
 787 eadd( ehalf, ex, ex );\r
 788 efloor( ex, ex ); /* y3 = floor( y3 ) - 27.0; */\r
 789 esub( e27, ex, en ); /* this is N */\r
 790 epow( eten, en, ex );\r
 791 emul( ex, em, erand );\r
 792 }\r
 793 \r
 794 /* -ln 2^16382 */\r
 795 char *ldemin = "-1.1355137111933024058873097E4";\r
 796 char *ldewid =  "2.2710274223866048117746193E4";\r
 797 /* -ln 2^1022 */\r
 798 char *demin  = "-7.0839641853226410622441123E2";\r
 799 char *dewid  =  "1.4167928370645282124488225E3";\r
 800 /* -ln 2^125 */\r
 801 char *femin  = "-8.6643397569993163677154015E1";\r
 802 char *fewid  =  "1.7328679513998632735430803E2";\r
 803 \r
 804 void etrand( erand )\r
 805 unsigned short erand[];\r
 806 {\r
 807 unsigned short ea[NE], ex[NE];\r
 808 double x, a;\r
 809 \r
 810 if( etrflag )\r
 811         {\r
 812         if( etrflag == 3 )\r
 813                 {\r
 814                 asctoe( ldemin, LOW );\r
 815                 asctoe( ldewid, WIDTH );\r
 816                 asctoe( a34, e27 );\r
 817                 }\r
 818         if( etrflag == 2 )\r
 819                 {\r
 820                 asctoe( demin, LOW );\r
 821                 asctoe( dewid, WIDTH );\r
 822                 asctoe( a27, e27 );\r
 823                 }\r
 824         if( etrflag == 1 )\r
 825                 {\r
 826                 asctoe( femin, LOW );\r
 827                 asctoe( fewid, WIDTH );\r
 828                 asctoe( a13, e27 );\r
 829                 }\r
 830         asctoe( a10m13, e10m13 );\r
 831         etrflag = 0;\r
 832         }\r
 833 drand( &x );\r
 834 e53toe( &x, ex ); /* x = WIDTH *  ( x - 1.0 )  +  LOW; */\r
 835 esub( eone, ex, ex );\r
 836 emul( WIDTH, ex, ex );\r
 837 eadd( LOW, ex, ex );\r
 838 eexp( ex, ex );   /* x = exp(x); */\r
 839 \r
 840 /* add fuzz\r
 841  */\r
 842 drand( &a );\r
 843 e53toe( &a, ea );\r
 844 emul( ea, ex, ea );  /* a = 1.0e-13 * x * a; */\r
 845 emul( e10m13, ea, ea );\r
 846 if( ecmp( ex, ezero ) > 0 )\r
 847         eneg( ea );\r
 848 eadd( ea, ex, erand );\r
 849 }\r
 850 \r