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
main()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 */
chkinternal(ref,tst,string)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 */
chkid(print,scan,string)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 */
chkscan(ref,tst,string)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 */
chkprint(ref,tst,string)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 */
pvec(x,n)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 */
cmpprint(ref,tst)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 */
cmpscan(ref,tst)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 */
cmptrunc(ref,tst)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
shownoncrit()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 */
chkit(testno)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
printerr()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
mnrand(erand)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
etrand(erand)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