Add in some math lib tests

1 files changed, 850 insertions, 0 deletions

diff --git a/test/math/ieetst.c b/test/math/ieetst.c
new file mode 100644
index 000000000..3c89145e3
--- /dev/null
+++ b/test/math/ieetst.c
@@ -0,0 +1,850 @@
+/* Floating point to ASCII input and output string test program.

+ *

+ * Numbers in the native machine data structure are converted

+ * to e type, then to and from decimal ASCII strings.  Native

+ * printf() and scanf() functions are also used to produce

+ * and read strings.  The resulting e type binary values

+ * are compared, with diagnostic printouts of any discrepancies.

+ *

+ * Steve Moshier, 16 Dec 88

+ * last revision: 16 May 92

+ */

+

+#include "ehead.h"

+#include "mconf.h"

+

+/* Include tests of 80-bit long double precision: */

+#define LDOUBLE 0

+/* Abort subtest after getting this many errors: */

+#define MAXERR 5

+/* Number of random arguments to try (set as large as you have

+ * patience for): */

+#define NRAND 100

+/* Perform internal consistency test: */

+#define CHKINTERNAL 0

+

+static unsigned short fullp[NE], rounded[NE];

+float prec24, sprec24, ssprec24;

+double prec53, sprec53, ssprec53;

+#if LDOUBLE

+long double prec64, sprec64, ssprec64;

+#endif

+

+static unsigned short rprint[NE], rscan[NE];

+static unsigned short q1[NE], q2[NE], q5[NE];

+static unsigned short e1[NE], e2[NE], e3[NE];

+static double d1, d2;

+static int errprint = 0;

+static int errscan = 0;

+static int identerr = 0;

+static int errtot = 0;

+static int count = 0;

+static char str0[80], str1[80], str2[80], str3[80];

+static unsigned short eten[NE], maxm[NE];

+

+int m, n, k2, mprec, SPREC;

+

+char *Ten = "10.0";

+char tformat[10];

+char *format24 = "%.8e";

+#ifdef DEC

+char *format53 = "%.17e";

+#else

+char *format53 = "%.16e";

+#endif

+char *fformat24 = "%e";

+char *fformat53 = "%le";

+char *pct = "%";

+char *quo = "\042";

+#if LDOUBLE

+char *format64 = "%.20Le";

+char *fformat64 = "%Le";

+#endif

+char *format;

+char *fformat;

+char *toomany = "Too many errors; aborting this test.\n";

+

+static int mnrflag;

+static int etrflag;

+void chkit(), printerr(), mnrand(), etrand(), shownoncrit();

+void chkid(), pvec();

+

+main()

+{

+int i, iprec;

+

+printf( "Steve Moshier's printf/scanf tester, version 0.2.\n\n" );

+#ifdef DEC

+ /* DEC PDP-11/VAX single precision not yet implemented */

+for( iprec = 1; iprec<2; iprec++ )

+#else

+for( iprec = 0; iprec<3; iprec++ )

+#endif

+	{

+	errscan = 0;

+	identerr = 0;

+	errprint = 0;

+	eclear( rprint );

+	eclear( rscan );

+

+switch( iprec )

+	{

+	case 0:

+		SPREC = 8; /* # digits after the decimal point */

+		mprec = 24; /* # bits in the significand */

+		m = 9; /* max # decimal digits for correct rounding */

+		n = 13; /* max power of ten for correct rounding */

+		k2 = -125; /* underflow beyond 2^-k2 */

+		format = format24; /* printf format string */

+		fformat = fformat24; /* scanf format string */

+		mnrflag = 1; /* sets interval for random numbers */

+		etrflag = 1;

+		printf( "Testing FLOAT precision.\n" );

+		break;

+

+	case 1:

+#ifdef DEC

+		SPREC = 17;

+		mprec = 56;

+		m = 17;

+		n = 27;

+		k2 = -125;

+		format = format53;

+		fformat = fformat53;

+		mnrflag = 2;

+		etrflag = 1;

+		printf( "Testing DEC DOUBLE precision.\n" );

+		break;

+#else

+		SPREC = 16;

+		mprec = 53;

+		m = 17;

+		n = 27;

+		k2 = -1021;

+		format = format53;

+		fformat = fformat53;

+		mnrflag = 2;

+		etrflag = 2;

+		printf( "Testing DOUBLE precision.\n" );

+		break;

+#endif

+	case 2:

+#if LDOUBLE

+		SPREC = 20;

+		mprec = 64;

+		m = 20;

+		n = 34;

+		k2 = -16382;

+		format = format64;

+		fformat = fformat64;

+		mnrflag = 3;

+		etrflag = 3;

+		printf( "Testing LONG DOUBLE precision.\n" );

+		break;

+#else

+		goto nodenorm;

+#endif

+	}

+

+	asctoe( Ten, eten );

+/* 10^m - 1 */

+	d2 = m;

+	e53toe( &d2, e1 );

+	epow( eten, e1, maxm );

+	esub( eone, maxm, maxm );

+

+/* test 1 */

+	printf( "1. Checking 10^n - 1 for n = %d to %d.\n", -m, m );

+	emov( eone, q5 );

+	for( count=0; count<=m; count++ )

+		{

+		esub( eone, q5, fullp );

+		chkit( 1 );

+		ediv( q5, eone, q2 );

+		esub( eone, q2, fullp );

+		chkit( 1 );

+		emul( eten, q5, q5 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end1;

+			}

+		}

+end1:

+	printerr();

+

+

+/* test 2 */

+	printf( "2. Checking powers of 10 from 10^-%d to 10^%d.\n", n, n );

+	emov( eone, q5 );

+	for( count=0; count<=n; count++ )

+		{

+		emov( q5, fullp );

+		chkit( 2 );

+		ediv( q5, eone, fullp );

+		chkit( 2 );

+		emul( eten, q5, q5 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end2;

+			}

+		}

+end2:

+	printerr();

+

+/* test 3 */

+	printf( "3. Checking (10^%d-1)*10^n from n = -%d to %d.\n", m, n, n );

+	emov( eone, q5 );

+	for( count= -n; count<=n; count++ )

+		{

+		emul( maxm, q5, fullp );

+		chkit( 3 );

+		emov( q5, fullp );

+		ediv( fullp, eone, fullp );

+		emul( maxm, fullp, fullp );

+		chkit( 3 );

+		emul( eten, q5, q5 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end3;

+			}

+		}

+end3:

+	printerr();

+

+

+

+/* test 4 */

+	printf( "4. Checking powers of 2 from 2^-24 to 2^+56.\n" );

+	d1 = -24.0;

+	e53toe( &d1, q1 );

+	epow( etwo, q1, q5 );

+

+	for( count = -24; count <= 56; count++ )

+		{

+		emov( q5, fullp );

+		chkit( 4 );

+		emul( etwo, q5, q5 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end4;

+			}

+		}

+end4:

+	printerr();

+

+

+/* test 5 */

+	printf( "5. Checking 2^n - 1 for n = 0 to %d.\n", mprec );

+	emov( eone, q5 );

+	for( count=0; count<=mprec; count++ )

+		{

+		esub( eone, q5, fullp );

+		chkit( 5 );

+		emul( etwo, q5, q5 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end5;

+			}

+		}

+end5:

+	printerr();

+

+/* test 6 */

+	printf( "6. Checking 2^n + 1 for n = 0 to %d.\n", mprec );

+	emov( eone, q5 );

+	for( count=0; count<=mprec; count++ )

+		{

+		eadd( eone, q5, fullp );

+		chkit( 6 );

+		emul( etwo, q5, q5 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end6;

+			}

+		}

+end6:

+	printerr();

+

+/* test 7 */

+	printf(

+	 "7. Checking %d values M * 10^N with random integer M and N,\n",

+	 NRAND );

+	printf("  1 <= M <= 10^%d - 1  and  -%d <= N <= +%d.\n", m, n, n );

+	for( i=0; i<NRAND; i++ )

+		{

+		mnrand( fullp );

+		chkit( 7 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end7;

+			}

+		}

+end7:

+	printerr();

+

+/* test 8 */

+	printf("8. Checking critical rounding cases.\n" );

+	for( i=0; i<20; i++ )

+		{

+		mnrand( fullp );

+		eabs( fullp );

+		if( ecmp( fullp, eone ) < 0 )

+			ediv( fullp, eone, fullp );

+		efloor( fullp, fullp );

+		eadd( ehalf, fullp, fullp );

+		chkit( 8 );

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto end8;

+			}

+		}

+end8:

+	printerr();

+

+

+

+/* test 9 */

+	printf("9. Testing on %d random non-denormal values.\n", NRAND );

+	for( i=0; i<NRAND; i++ )

+		{

+		etrand( fullp );

+		chkit( 9 );

+		}

+	printerr();

+	shownoncrit();

+

+/* test 10 */

+	printf(

+	"Do you want to check denormal numbers in this precision ? (y/n) " );

+	gets( str0 );

+	if( str0[0] != 'y' )

+		goto nodenorm;

+

+	printf( "10. Checking denormal numbers.\n" );

+

+/* Form 2^-starting power */

+	d1 = k2;

+	e53toe( &d1, q1 );

+	epow( etwo, q1, e1 );

+

+/* Find 2^-mprec less than starting power */

+	d1 = -mprec + 4;

+	e53toe( &d1, q1 );

+	epow( etwo, q1, e3 );

+	emul( e1, e3, e3 );

+	emov( e3, e2 );

+	ediv( etwo, e2, e2 );

+

+	while( ecmp(e1,e2) != 0 )

+		{

+		eadd( e1, e2, fullp );

+		switch( mprec )

+			{

+#if LDOUBLE

+			case 64:

+			etoe64( e1, &sprec64 );

+			e64toe( &sprec64, q1 );

+			etoe64( fullp, &prec64 );

+			e64toe( &prec64, q2 );

+			break;

+#endif

+#ifdef DEC

+			case 56:

+#endif

+			case 53:

+			etoe53( e1, &sprec53 );

+			e53toe( &sprec53, q1 );

+			etoe53( fullp, &prec53 );

+			e53toe( &prec53, q2 );

+			break;

+

+			case 24:

+			etoe24( e1, &sprec24 );

+			e24toe( &sprec24, q1 );

+			etoe24( fullp, &prec24 );

+			e24toe( &prec24, q2 );

+			break;

+			}

+		if( ecmp( q2, ezero ) == 0 )

+			goto maxden;

+		chkit(10);

+		if( ecmp(q1,q2) == 0 )

+			{

+			ediv( etwo, e1, e1 );

+			emov( e3, e2 );

+			}

+		if( errtot >= MAXERR )

+			{

+			printf( "%s", toomany );

+			goto maxden;

+			}

+		ediv( etwo, e2, e2 );

+		}

+maxden:

+	printerr();

+nodenorm:

+	printf( "\n" );

+	} /* loop on precision */

+printf( "End of test.\n" );

+}

+

+#if CHKINTERNAL

+long double xprec64;

+double xprec53;

+float xprec24;

+

+/* Check binary -> printf -> scanf -> binary identity

+ * of internal routines

+ */

+void chkinternal( ref, tst, string )

+unsigned short ref[], tst[];

+char *string;

+{

+

+if( ecmp(ref,tst) != 0 )

+	{

+	printf( "internal identity compare error!\n" );

+	chkid( ref, tst, string );

+	}

+}

+#endif

+

+

+/* Check binary -> printf -> scanf -> binary identity

+ */

+void chkid( print, scan, string )

+unsigned short print[], scan[];

+char *string;

+{

+/* Test printf-scanf identity */

+if( ecmp( print, scan ) != 0 )

+	{

+	pvec( print, NE );

+	printf( " ->printf-> %s ->scanf->\n", string );

+	pvec( scan, NE );

+	printf( " is not an identity.\n" );

+	++identerr;

+	}

+}

+

+

+/* Check scanf result

+ */

+void chkscan( ref, tst, string )

+unsigned short ref[], tst[];

+char *string;

+{

+/* Test scanf()  */

+if( ecmp( ref, tst ) != 0 )

+	{

+	printf( "scanf(%s) -> ", string );

+	pvec( tst, NE );

+	printf( "\n should be    " );

+	pvec( ref, NE );

+	printf( ".\n" );

+	++errscan;

+	++errtot;

+	}

+}

+

+

+/* Test printf() result

+ */

+void chkprint( ref, tst, string ) 

+unsigned short ref[], tst[];

+char *string;

+{

+if( ecmp(ref, tst) != 0 )

+	{

+	printf( "printf( ");

+	pvec( ref, NE );

+	printf( ") -> %s\n", string );

+	printf( "      = " );

+	pvec( tst, NE );

+	printf( ".\n" );

+	++errprint;

+	++errtot;

+	}

+}

+

+

+/* Print array of n 16-bit shorts

+ */

+void pvec( x, n )

+unsigned short x[];

+int n;

+{

+int i;

+

+for( i=0; i<n; i++ )

+	{

+	printf( "%04x ", x[i] );

+	}

+}

+

+/* Measure worst case printf rounding error

+ */

+void cmpprint( ref, tst )

+unsigned short ref[], tst[];

+{

+unsigned short e[NE];

+

+if( ecmp( ref, ezero ) != 0 )

+	{

+	esub( ref, tst, e );

+	ediv( ref, e, e );

+	eabs( e );

+	if( ecmp( e, rprint ) > 0 )

+		emov( e, rprint );

+	}

+}

+

+/* Measure worst case scanf rounding error

+ */

+void cmpscan( ref, tst )

+unsigned short ref[], tst[];

+{

+unsigned short er[NE];

+

+if( ecmp( ref, ezero ) != 0 )

+	{

+	esub( ref, tst, er );

+	ediv( ref, er, er );

+	eabs( er );

+	if( ecmp( er, rscan ) > 0 )

+		emov( er, rscan );

+	if( ecmp( er, ehalf ) > 0 )

+		{

+		etoasc( tst, str1, 21 );

+		printf( "Bad error: scanf(%s) = %s !\n", str0, str1 );

+		}

+	}

+}

+

+/* Check rounded-down decimal string output of printf

+ */

+void cmptrunc( ref, tst )

+unsigned short ref[], tst[];

+{

+if( ecmp( ref, tst ) != 0 )

+	{

+	printf( "printf(%s%s%s, %s) -> %s\n", quo, tformat, quo, str1, str2 );

+	printf( "should be      %s .\n", str3 );

+	errprint += 1;

+	}

+}

+

+

+void shownoncrit()

+{

+

+etoasc( rprint, str0, 3 );

+printf( "Maximum relative printf error found = %s .\n", str0 );

+etoasc( rscan, str0, 3 );

+printf( "Maximum relative scanf error found = %s .\n", str0 );

+}

+

+

+

+/* Produce arguments and call comparison subroutines.

+ */

+void chkit( testno )

+int testno;

+{

+unsigned short t[NE], u[NE], v[NE];

+int j;

+

+switch( mprec )

+	{

+#if LDOUBLE

+	case 64:

+		etoe64( fullp, &prec64 );

+		e64toe( &prec64, rounded );

+#if CHKINTERNAL

+		e64toasc( &prec64, str1, SPREC );

+		asctoe64( str1, &xprec64 );

+		e64toe( &xprec64, t );

+		chkinternal( rounded, t, str1 );

+#endif

+/* check printf and scanf */

+		sprintf( str2, format, prec64 );

+		sscanf( str2, fformat, &sprec64 );

+		e64toe( &sprec64, u );

+		chkid( rounded, u, str2 );

+		asctoe64( str2, &ssprec64 );

+		e64toe( &ssprec64, v );

+		chkscan( v, u, str2 );

+		chkprint( rounded, v, str2 );

+		if( testno < 8 )

+			break;

+/* rounding error measurement */

+		etoasc( fullp, str0, 24 );

+		etoe64( fullp, &ssprec64 );

+		e64toe( &ssprec64, u );

+		sprintf( str2, format, ssprec64 );

+		asctoe( str2, t );

+		cmpprint( u, t );

+		sscanf( str0, fformat, &sprec64 );

+		e64toe( &sprec64, t );

+		cmpscan( fullp, t );

+		if( testno < 8 )

+			break;

+/* strings rounded to less than maximum precision */

+		e64toasc( &ssprec64, str1, 24 );

+		for( j=SPREC-1; j>0; j-- )		

+			{

+			e64toasc( &ssprec64, str3, j );

+			asctoe( str3, v );

+			sprintf( tformat, "%s.%dLe", pct, j );

+			sprintf( str2, tformat, ssprec64 );

+			asctoe( str2, t );

+			cmptrunc( v, t );

+			}

+		break;

+#endif

+#ifdef DEC

+	case 56:

+#endif

+	case 53:

+		etoe53( fullp, &prec53 );

+		e53toe( &prec53, rounded );

+#if CHKINTERNAL

+		e53toasc( &prec53, str1, SPREC );

+		asctoe53( str1, &xprec53 );

+		e53toe( &xprec53, t );

+		chkinternal( rounded, t, str1 );

+#endif

+		sprintf( str2, format, prec53 );

+		sscanf( str2, fformat, &sprec53 );

+		e53toe( &sprec53, u );

+		chkid( rounded, u, str2 );

+		asctoe53( str2, &ssprec53 );

+		e53toe( &ssprec53, v );

+		chkscan( v, u, str2 );

+		chkprint( rounded, v, str2 );

+		if( testno < 8 )

+			break;

+/* rounding error measurement */

+		etoasc( fullp, str0, 24 );

+		etoe53( fullp, &ssprec53 );

+		e53toe( &ssprec53, u );

+		sprintf( str2, format, ssprec53 );

+		asctoe( str2, t );

+		cmpprint( u, t );

+		sscanf( str0, fformat, &sprec53 );

+		e53toe( &sprec53, t );

+		cmpscan( fullp, t );

+		if( testno < 8 )

+			break;

+		e53toasc( &ssprec53, str1, 24 );

+		for( j=SPREC-1; j>0; j-- )		

+			{

+			e53toasc( &ssprec53, str3, j );

+			asctoe( str3, v );

+			sprintf( tformat, "%s.%de", pct, j );

+			sprintf( str2, tformat, ssprec53 );

+			asctoe( str2, t );

+			cmptrunc( v, t );

+			}

+		break;

+

+	case 24:

+		etoe24( fullp, &prec24 );

+		e24toe( &prec24, rounded );

+#if CHKINTERNAL

+		e24toasc( &prec24, str1, SPREC );

+		asctoe24( str1, &xprec24 );

+		e24toe( &xprec24, t );

+		chkinternal( rounded, t, str1 );

+#endif

+		sprintf( str2, format, prec24 );

+		sscanf( str2, fformat, &sprec24 );

+		e24toe( &sprec24, u );

+		chkid( rounded, u, str2 );

+		asctoe24( str2, &ssprec24 );

+		e24toe( &ssprec24, v );

+		chkscan( v, u, str2 );

+		chkprint( rounded, v, str2 );

+		if( testno < 8 )

+			break;

+/* rounding error measurement */

+		etoasc( fullp, str0, 24 );

+		etoe24( fullp, &ssprec24 );

+		e24toe( &ssprec24, u );

+		sprintf( str2, format, ssprec24 );

+		asctoe( str2, t );

+		cmpprint( u, t );

+		sscanf( str0, fformat, &sprec24 );

+		e24toe( &sprec24, t );

+		cmpscan( fullp, t );

+/*

+		if( testno < 8 )

+			break;

+*/

+		e24toasc( &ssprec24, str1, 24 );

+		for( j=SPREC-1; j>0; j-- )		

+			{

+			e24toasc( &ssprec24, str3, j );

+			asctoe( str3, v );

+			sprintf( tformat, "%s.%de", pct, j );

+			sprintf( str2, tformat, ssprec24 );

+			asctoe( str2, t );

+			cmptrunc( v, t );

+			}

+		break;

+	}

+}

+

+

+void printerr()

+{

+if( (errscan == 0) && (identerr == 0) && (errprint == 0) )

+	printf( "No errors found.\n" );

+else

+	{

+	printf( "%d binary -> decimal errors found.\n", errprint );

+	printf( "%d decimal -> binary errors found.\n", errscan );

+	}

+errscan = 0;	/* reset for next test */

+identerr = 0;

+errprint = 0;

+errtot = 0;

+}

+

+

+/* Random number generator

+ * in the range M * 10^N, where 1 <= M <= 10^17 - 1

+ * and -27 <= N <= +27.  Test values of M are logarithmically distributed

+ * random integers; test values of N are uniformly distributed random integers.

+ */

+

+static char *fwidth = "1.036163291797320557783096e1"; /* log(sqrt(10^9-1)) */

+static char *dwidth = "1.957197329044938830915E1"; /* log(sqrt(10^17-1)) */

+static char *ldwidth = "2.302585092994045684017491e1"; /* log(sqrt(10^20-1)) */

+

+static char *a13 = "13.0";

+static char *a27 = "27.0";

+static char *a34 = "34.0";

+static char *a10m13 = "1.0e-13";

+static unsigned short LOW[ NE ], WIDTH[NE], e27[NE], e10m13[NE];

+

+

+void mnrand( erand )

+unsigned short erand[];

+{

+unsigned short ea[NE], em[NE], en[NE], ex[NE];

+double x, a;

+

+if( mnrflag )

+	{

+	if( mnrflag == 3 )

+		{

+		asctoe( ldwidth, WIDTH );

+		asctoe( a34, e27 );

+		}

+	if( mnrflag == 2 )

+		{

+		asctoe( dwidth, WIDTH );

+		asctoe( a27, e27 );

+		}

+	if( mnrflag == 1 )

+		{

+		asctoe( fwidth, WIDTH );

+		asctoe( a13, e27 );

+		}

+	asctoe( a10m13, e10m13 );

+	mnrflag = 0;

+	}

+drand( &x );

+e53toe( &x, ex ); /* x = WIDTH *  ( x - 1.0 )  +  LOW; */

+esub( eone, ex, ex );

+emul( WIDTH, ex, ex );

+eexp( ex, ex );   /* x = exp(x); */

+

+drand( &a );

+e53toe( &a, ea );

+emul( ea, ex, ea );  /* a = 1.0e-13 * x * a; */

+emul( e10m13, ea, ea );

+eabs( ea );

+eadd( ea, ex, ex );	/* add fuzz */

+emul( ex, ex, ex );	/* square it, to get range to 10^17 - 1 */

+efloor( ex, em ); /* this is M */

+

+/* Random power of 10 */

+drand( &a );

+e53toe( &a, ex );

+esub( eone, ex, ex ); /* y3 = 54.0 *  ( y3 - 1.0 ) + 0.5; */

+emul( e27, ex, ex );

+eadd( ex, ex, ex );

+eadd( ehalf, ex, ex );

+efloor( ex, ex ); /* y3 = floor( y3 ) - 27.0; */

+esub( e27, ex, en ); /* this is N */

+epow( eten, en, ex );

+emul( ex, em, erand );

+}

+

+/* -ln 2^16382 */

+char *ldemin = "-1.1355137111933024058873097E4";

+char *ldewid =  "2.2710274223866048117746193E4";

+/* -ln 2^1022 */

+char *demin  = "-7.0839641853226410622441123E2";

+char *dewid  =  "1.4167928370645282124488225E3";

+/* -ln 2^125 */

+char *femin  = "-8.6643397569993163677154015E1";

+char *fewid  =  "1.7328679513998632735430803E2";

+

+void etrand( erand )

+unsigned short erand[];

+{

+unsigned short ea[NE], ex[NE];

+double x, a;

+

+if( etrflag )

+	{

+	if( etrflag == 3 )

+		{

+		asctoe( ldemin, LOW );

+		asctoe( ldewid, WIDTH );

+		asctoe( a34, e27 );

+		}

+	if( etrflag == 2 )

+		{

+		asctoe( demin, LOW );

+		asctoe( dewid, WIDTH );

+		asctoe( a27, e27 );

+		}

+	if( etrflag == 1 )

+		{

+		asctoe( femin, LOW );

+		asctoe( fewid, WIDTH );

+		asctoe( a13, e27 );

+		}

+	asctoe( a10m13, e10m13 );

+	etrflag = 0;

+	}

+drand( &x );

+e53toe( &x, ex ); /* x = WIDTH *  ( x - 1.0 )  +  LOW; */

+esub( eone, ex, ex );

+emul( WIDTH, ex, ex );

+eadd( LOW, ex, ex );

+eexp( ex, ex );   /* x = exp(x); */

+

+/* add fuzz

+ */

+drand( &a );

+e53toe( &a, ea );

+emul( ea, ex, ea );  /* a = 1.0e-13 * x * a; */

+emul( e10m13, ea, ea );

+if( ecmp( ex, ezero ) > 0 )

+	eneg( ea );

+eadd( ea, ex, erand );

+}

+