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