summaryrefslogtreecommitdiff
path: root/src/regress/lib/libc/cephes/ieetst.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/regress/lib/libc/cephes/ieetst.c')
-rw-r--r--src/regress/lib/libc/cephes/ieetst.c880
1 files changed, 0 insertions, 880 deletions
diff --git a/src/regress/lib/libc/cephes/ieetst.c b/src/regress/lib/libc/cephes/ieetst.c
deleted file mode 100644
index 41dfd6a96e..0000000000
--- a/src/regress/lib/libc/cephes/ieetst.c
+++ /dev/null
@@ -1,880 +0,0 @@
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
51static unsigned short fullp[NE], rounded[NE];
52float prec24, sprec24, ssprec24;
53double prec53, sprec53, ssprec53;
54#if LDOUBLE
55long double prec64, sprec64, ssprec64;
56#endif
57
58static unsigned short rprint[NE], rscan[NE];
59static unsigned short q1[NE], q2[NE], q5[NE];
60static unsigned short e1[NE], e2[NE], e3[NE];
61static double d1, d2;
62static int errprint = 0;
63static int errscan = 0;
64static int identerr = 0;
65static int errtot = 0;
66static int count = 0;
67static char str0[80], str1[80], str2[80], str3[80];
68static unsigned short eten[NE], maxm[NE];
69
70int m, n, k2, mprec, SPREC;
71
72char *Ten = "10.0";
73char tformat[10];
74char *format24 = "%.8e";
75#ifdef DEC
76char *format53 = "%.17e";
77#else
78char *format53 = "%.16e";
79#endif
80char *fformat24 = "%e";
81char *fformat53 = "%le";
82char *pct = "%";
83char *quo = "\042";
84#if LDOUBLE
85char *format64 = "%.20Le";
86char *fformat64 = "%Le";
87#endif
88char *format;
89char *fformat;
90char *toomany = "Too many errors; aborting this test.\n";
91
92static int mnrflag;
93static int etrflag;
94void chkit(), printerr(), mnrand(), etrand(), shownoncrit();
95void chkid(), pvec();
96
97main()
98{
99int i, iprec, retval = 0;
100
101printf( "Steve Moshier's printf/scanf tester, version 0.2.\n\n" );
102#ifdef DEC
103 /* DEC PDP-11/VAX single precision not yet implemented */
104for( iprec = 1; iprec<2; iprec++ )
105#else
106for( 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
116switch( 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 }
199end1:
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 }
219end2:
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 }
240end3:
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 }
262end4:
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 }
280end5:
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 }
297end6:
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 }
315end7:
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 }
335end8:
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 }
419maxden:
420 printerr();
421nodenorm:
422 printf( "\n" );
423 retval |= errscan | identerr | errprint;
424 } /* loop on precision */
425printf( "End of test.\n" );
426return (retval);
427}
428
429#if CHKINTERNAL
430long double xprec64;
431double xprec53;
432float xprec24;
433
434/* Check binary -> printf -> scanf -> binary identity
435 * of internal routines
436 */
437void chkinternal( ref, tst, string )
438unsigned short ref[], tst[];
439char *string;
440{
441
442if( 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 */
453void chkid( print, scan, string )
454unsigned short print[], scan[];
455char *string;
456{
457/* Test printf-scanf identity */
458if( 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 */
471void chkscan( ref, tst, string )
472unsigned short ref[], tst[];
473char *string;
474{
475/* Test scanf() */
476if( 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 */
491void chkprint( ref, tst, string )
492unsigned short ref[], tst[];
493char *string;
494{
495if( 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 */
511void pvec( x, n )
512unsigned short x[];
513int n;
514{
515int i;
516
517for( i=0; i<n; i++ )
518 {
519 printf( "%04x ", x[i] );
520 }
521}
522
523/* Measure worst case printf rounding error
524 */
525void cmpprint( ref, tst )
526unsigned short ref[], tst[];
527{
528unsigned short e[NE];
529
530if( 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 */
542void cmpscan( ref, tst )
543unsigned short ref[], tst[];
544{
545unsigned short er[NE];
546
547if( 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 */
564void cmptrunc( ref, tst )
565unsigned short ref[], tst[];
566{
567if( 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
576void shownoncrit()
577{
578
579etoasc( rprint, str0, 3 );
580printf( "Maximum relative printf error found = %s .\n", str0 );
581etoasc( rscan, str0, 3 );
582printf( "Maximum relative scanf error found = %s .\n", str0 );
583}
584
585
586
587/* Produce arguments and call comparison subroutines.
588 */
589void chkit( testno )
590int testno;
591{
592unsigned short t[NE], u[NE], v[NE];
593int j;
594
595switch( 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
737void printerr()
738{
739if( (errscan == 0) && (identerr == 0) && (errprint == 0) )
740 printf( "No errors found.\n" );
741else
742 {
743 printf( "%d binary -> decimal errors found.\n", errprint );
744 printf( "%d decimal -> binary errors found.\n", errscan );
745 }
746errscan = 0; /* reset for next test */
747identerr = 0;
748errprint = 0;
749errtot = 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
759static char *fwidth = "1.036163291797320557783096e1"; /* log(sqrt(10^9-1)) */
760static char *dwidth = "1.957197329044938830915E1"; /* log(sqrt(10^17-1)) */
761static char *ldwidth = "2.302585092994045684017491e1"; /* log(sqrt(10^20-1)) */
762
763static char *a13 = "13.0";
764static char *a27 = "27.0";
765static char *a34 = "34.0";
766static char *a10m13 = "1.0e-13";
767static unsigned short LOW[ NE ], WIDTH[NE], e27[NE], e10m13[NE];
768
769
770void mnrand( erand )
771unsigned short erand[];
772{
773unsigned short ea[NE], em[NE], en[NE], ex[NE];
774double x, a;
775
776if( 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 }
796drand( &x );
797e53toe( &x, ex ); /* x = WIDTH * ( x - 1.0 ) + LOW; */
798esub( eone, ex, ex );
799emul( WIDTH, ex, ex );
800eexp( ex, ex ); /* x = exp(x); */
801
802drand( &a );
803e53toe( &a, ea );
804emul( ea, ex, ea ); /* a = 1.0e-13 * x * a; */
805emul( e10m13, ea, ea );
806eabs( ea );
807eadd( ea, ex, ex ); /* add fuzz */
808emul( ex, ex, ex ); /* square it, to get range to 10^17 - 1 */
809efloor( ex, em ); /* this is M */
810
811/* Random power of 10 */
812drand( &a );
813e53toe( &a, ex );
814esub( eone, ex, ex ); /* y3 = 54.0 * ( y3 - 1.0 ) + 0.5; */
815emul( e27, ex, ex );
816eadd( ex, ex, ex );
817eadd( ehalf, ex, ex );
818efloor( ex, ex ); /* y3 = floor( y3 ) - 27.0; */
819esub( e27, ex, en ); /* this is N */
820epow( eten, en, ex );
821emul( ex, em, erand );
822}
823
824/* -ln 2^16382 */
825char *ldemin = "-1.1355137111933024058873097E4";
826char *ldewid = "2.2710274223866048117746193E4";
827/* -ln 2^1022 */
828char *demin = "-7.0839641853226410622441123E2";
829char *dewid = "1.4167928370645282124488225E3";
830/* -ln 2^125 */
831char *femin = "-8.6643397569993163677154015E1";
832char *fewid = "1.7328679513998632735430803E2";
833
834void etrand( erand )
835unsigned short erand[];
836{
837unsigned short ea[NE], ex[NE];
838double x, a;
839
840if( 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 }
863drand( &x );
864e53toe( &x, ex ); /* x = WIDTH * ( x - 1.0 ) + LOW; */
865esub( eone, ex, ex );
866emul( WIDTH, ex, ex );
867eadd( LOW, ex, ex );
868eexp( ex, ex ); /* x = exp(x); */
869
870/* add fuzz
871 */
872drand( &a );
873e53toe( &a, ea );
874emul( ea, ex, ea ); /* a = 1.0e-13 * x * a; */
875emul( e10m13, ea, ea );
876if( ecmp( ex, ezero ) > 0 )
877 eneg( ea );
878eadd( ea, ex, erand );
879}
880
generated by cgit v1.2.3-55-g6feb (git 2.39.1) at 2025-11-19 23:39:19 +0000