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