Add in some math lib tests

15 files changed, 9839 insertions, 0 deletions

diff --git a/test/math/Makefile b/test/math/Makefile
new file mode 100644
index 000000000..f8321b92e
--- /dev/null
+++ b/test/math/Makefile
@@ -0,0 +1,71 @@
+# Unix makefile for ieetst, eparanoi.

+# Set LARGEMEM 1 in qcalc.h for 32-bit memory addresses.

+# Define computer type and/or endianness in mconf.h.

+#

+# Configure eparanoi.c for desired arithmetic test;

+# also define appropriate version of setprec.o, or use a stub that

+# does no FPU setup.  To test native arithmetic, eparanoi uses

+# the system libraries only; compile simply by `cc eparanoi.c -lm'.

+#

+

+TESTDIR=../

+include $(TESTDIR)/Rules.mak

+

+

+#CC = gcc

+#CFLAGS= -O

+INCS= mconf.h ehead.h

+OBJS = ieee.o econst.o eexp.o elog.o epow.o etanh.o etodec.o mtherr.o #setprec.o

+TARGETS=ieetst eparanoi

+

+all: $(TARGETS)

+

+ieetst: ieetst.o $(OBJS) drand.o $(INCS)

+	$(CC) -o ieetst ieetst.o $(OBJS) drand.o -lc -lm

+

+eparanoi: eparanoi.o $(OBJS) $(INCS)

+	$(CC) -o eparanoi  eparanoi.o $(OBJS) -lc -lm

+

+#setprec.o: setprec.387

+#	as -o setprec.o setprec.387

+

+#setprec.o: setprec.688

+#	as -o setprec.o setprec.688

+

+ieee.o: ieee.c $(INCS)

+	$(CC) $(CFLAGS) -c ieee.c

+

+econst.o: econst.c $(INCS)

+	$(CC) $(CFLAGS) -c econst.c

+

+elog.o: elog.c $(INCS)

+	$(CC) $(CFLAGS) -c elog.c

+

+eexp.o: eexp.c $(INCS)

+	$(CC) $(CFLAGS) -c eexp.c

+

+etanh.o: etanh.c $(INCS)

+	$(CC) $(CFLAGS) -c etanh.c

+

+epow.o: epow.c $(INCS)

+	$(CC) $(CFLAGS) -c epow.c

+

+mtherr.o: mtherr.c $(INCS)

+	$(CC) $(CFLAGS) -c mtherr.c

+

+ieetst.o: ieetst.c $(INCS)

+	$(CC) $(CFLAGS) -c ieetst.c

+

+drand.o: drand.c $(INCS)

+	$(CC) $(CFLAGS) -c drand.c

+

+etodec.o: etodec.c $(INCS)

+	$(CC) $(CFLAGS) -c etodec.c

+

+eparanoi.o: eparanoi.c $(INCS)

+	$(CC) $(CFLAGS) -c eparanoi.c

+

+clean:

+	rm -f *.[oa] *~ core $(TARGETS)

+

+

diff --git a/test/math/drand.c b/test/math/drand.c
new file mode 100644
index 000000000..9eedf71fc
--- /dev/null
+++ b/test/math/drand.c
@@ -0,0 +1,158 @@
+/*							drand.c
+ *
+ *	Pseudorandom number generator
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double y, drand();
+ *
+ * drand( &y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Yields a random number 1.0 <= y < 2.0.
+ *
+ * The three-generator congruential algorithm by Brian
+ * Wichmann and David Hill (BYTE magazine, March, 1987,
+ * pp 127-8) is used. The period, given by them, is
+ * 6953607871644.
+ *
+ * Versions invoked by the different arithmetic compile
+ * time options DEC, IBMPC, and MIEEE, produce
+ * approximately the same sequences, differing only in the
+ * least significant bits of the numbers. The UNK option
+ * implements the algorithm as recommended in the BYTE
+ * article.  It may be used on all computers. However,
+ * the low order bits of a double precision number may
+ * not be adequately random, and may vary due to arithmetic
+ * implementation details on different computers.
+ *
+ * The other compile options generate an additional random
+ * integer that overwrites the low order bits of the double
+ * precision number.  This reduces the period by a factor of
+ * two but tends to overcome the problems mentioned.
+ *
+ */
+
+
+
+#include "mconf.h"
+
+
+/*  Three-generator random number algorithm
+ * of Brian Wichmann and David Hill
+ * BYTE magazine, March, 1987 pp 127-8
+ *
+ * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
+ */
+
+static int sx = 1;
+static int sy = 10000;
+static int sz = 3000;
+
+static union {
+ double d;
+ unsigned short s[4];
+} unkans;
+
+/* This function implements the three
+ * congruential generators.
+ */
+ 
+int ranwh()
+{
+int r, s;
+
+/*  sx = sx * 171 mod 30269 */
+r = sx/177;
+s = sx - 177 * r;
+sx = 171 * s - 2 * r;
+if( sx < 0 )
+	sx += 30269;
+
+
+/* sy = sy * 172 mod 30307 */
+r = sy/176;
+s = sy - 176 * r;
+sy = 172 * s - 35 * r;
+if( sy < 0 )
+	sy += 30307;
+
+/* sz = 170 * sz mod 30323 */
+r = sz/178;
+s = sz - 178 * r;
+sz = 170 * s - 63 * r;
+if( sz < 0 )
+	sz += 30323;
+/* The results are in static sx, sy, sz. */
+return 0;
+}
+
+/*	drand.c
+ *
+ * Random double precision floating point number between 1 and 2.
+ *
+ * C callable:
+ *	drand( &x );
+ */
+
+int drand( a )
+double *a;
+{
+unsigned short r;
+#ifdef DEC
+unsigned short s, t;
+#endif
+
+/* This algorithm of Wichmann and Hill computes a floating point
+ * result:
+ */
+ranwh();
+unkans.d = sx/30269.0  +  sy/30307.0  +  sz/30323.0;
+r = unkans.d;
+unkans.d -= r;
+unkans.d += 1.0;
+
+/* if UNK option, do nothing further.
+ * Otherwise, make a random 16 bit integer
+ * to overwrite the least significant word
+ * of unkans.
+ */
+#ifdef UNK
+/* do nothing */
+#else
+ranwh();
+r = sx * sy + sz;
+#endif
+
+#ifdef DEC
+/* To make the numbers as similar as possible
+ * in all arithmetics, the random integer has
+ * to be inserted 3 bits higher up in a DEC number.
+ * An alternative would be put it 3 bits lower down
+ * in all the other number types.
+ */
+s = unkans.s[2];
+t = s & 07;	/* save these bits to put in at the bottom */
+s &= 0177770;
+s |= (r >> 13) & 07;
+unkans.s[2] = s;
+t |= r << 3;
+unkans.s[3] = t;
+#endif
+
+#ifdef IBMPC
+unkans.s[0] = r;
+#endif
+
+#ifdef MIEEE
+unkans.s[3] = r;
+#endif
+
+*a = unkans.d;
+return 0;
+}
diff --git a/test/math/econst.c b/test/math/econst.c
new file mode 100644
index 000000000..21523e568
--- /dev/null
+++ b/test/math/econst.c
@@ -0,0 +1,96 @@
+/*							econst.c	*/

+/*  e type constants used by high precision check routines */

+

+#include "ehead.h"

+

+

+#if NE == 10

+/* 0.0 */

+unsigned short ezero[NE] =

+ {0x0000, 0x0000, 0x0000, 0x0000,

+  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};

+

+/* 5.0E-1 */

+unsigned short ehalf[NE] =

+ {0x0000, 0x0000, 0x0000, 0x0000,

+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};

+

+/* 1.0E0 */

+unsigned short eone[NE] =

+ {0x0000, 0x0000, 0x0000, 0x0000,

+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};

+

+/* 2.0E0 */

+unsigned short etwo[NE] =

+ {0x0000, 0x0000, 0x0000, 0x0000,

+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};

+

+/* 3.2E1 */

+unsigned short e32[NE] =

+ {0x0000, 0x0000, 0x0000, 0x0000,

+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};

+

+/* 6.93147180559945309417232121458176568075500134360255E-1 */

+unsigned short elog2[NE] =

+ {0x40f3, 0xf6af, 0x03f2, 0xb398,

+  0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};

+

+/* 1.41421356237309504880168872420969807856967187537695E0 */

+unsigned short esqrt2[NE] =

+ {0x1d6f, 0xbe9f, 0x754a, 0x89b3,

+  0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};

+

+/* 3.14159265358979323846264338327950288419716939937511E0 */

+unsigned short epi[NE] =

+ {0x2902, 0x1cd1, 0x80dc, 0x628b,

+  0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};

+  

+/* 5.7721566490153286060651209008240243104215933593992E-1 */

+unsigned short eeul[NE] = {

+0xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,};

+

+#else

+

+/* 0.0 */

+unsigned short ezero[NE] = {

+0, 0000000,0000000,0000000,0000000,0000000,};

+/* 5.0E-1 */

+unsigned short ehalf[NE] = {

+0, 0000000,0000000,0000000,0100000,0x3ffe,};

+/* 1.0E0 */

+unsigned short eone[NE] = {

+0, 0000000,0000000,0000000,0100000,0x3fff,};

+/* 2.0E0 */

+unsigned short etwo[NE] = {

+0, 0000000,0000000,0000000,0100000,0040000,};

+/* 3.2E1 */

+unsigned short e32[NE] = {

+0, 0000000,0000000,0000000,0100000,0040004,};

+/* 6.93147180559945309417232121458176568075500134360255E-1 */

+unsigned short elog2[NE] = {

+0xc9e4,0x79ab,0150717,0013767,0130562,0x3ffe,};

+/* 1.41421356237309504880168872420969807856967187537695E0 */

+unsigned short esqrt2[NE] = {

+0x597e,0x6484,0174736,0171463,0132404,0x3fff,};

+/* 2/sqrt(PI) =

+ * 1.12837916709551257389615890312154517168810125865800E0 */

+unsigned short eoneopi[NE] = {

+0x71d5,0x688d,0012333,0135202,0110156,0x3fff,};

+/* 3.14159265358979323846264338327950288419716939937511E0 */

+unsigned short epi[NE] = {

+0xc4c6,0xc234,0020550,0155242,0144417,0040000,};

+/* 5.7721566490153286060651209008240243104215933593992E-1 */

+unsigned short eeul[NE] = {

+0xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,};

+#endif

+extern unsigned short ezero[];

+extern unsigned short ehalf[];

+extern unsigned short eone[];

+extern unsigned short etwo[];

+extern unsigned short e32[];

+extern unsigned short elog2[];

+extern unsigned short esqrt2[];

+extern unsigned short eoneopi[];

+extern unsigned short epi[];

+extern unsigned short eeul[];

+

diff --git a/test/math/eexp.c b/test/math/eexp.c
new file mode 100644
index 000000000..dd2a64e1c
--- /dev/null
+++ b/test/math/eexp.c
@@ -0,0 +1,77 @@
+/*							xexp.c		*/

+/* exponential function check routine */

+/* by Stephen L. Moshier. */

+

+

+#include "ehead.h"

+

+/*

+extern int powinited;

+extern short maxposint[], maxnegint[];

+*/

+

+void eexp( x, y )

+unsigned short *x, *y;

+{

+unsigned short num[NE], den[NE], x2[NE];

+long i;

+unsigned short sign, expchk;

+

+/* range reduction theory: x = i + f, 0<=f<1;

+ * e**x = e**i * e**f 

+ * e**i = 2**(i/log 2).

+ * Let i/log2 = i1 + f1, 0<=f1<1.

+ * Then e**i = 2**i1 * 2**f1, so

+ * e**x = 2**i1 * e**(log 2 * f1) * e**f.

+ */

+/*

+if( powinited == 0 )

+	initpow();

+*/

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

+	{

+	emov( eone, y );

+	return;

+	}

+emov(x, x2);

+expchk = x2[NE-1];

+sign = expchk & 0x8000;

+x2[NE-1] &= 0x7fff;

+

+/* Test for excessively large argument */

+expchk &= 0x7fff;

+if( expchk > (EXONE + 15) )

+	{

+	eclear( y );

+	if( sign == 0 )

+		einfin( y );

+	return;

+	}

+

+eifrac( x2, &i, num );		/* x = i + f		*/

+

+if( i != 0 )

+ {

+ ltoe( &i, den );		/* floating point i	*/

+ ediv( elog2, den, den );	/* i/log 2		*/

+ eifrac( den, &i, den );	/* i/log 2  =  i1 + f1	*/

+ emul( elog2, den, den );	/* log 2 * f1		*/

+ eadd( den, num, x2 );		/* log 2 * f1  + f	*/

+ }

+

+/*x2[NE-1] -= 1;*/

+eldexp( x2, -1L, x2 ); /* divide by 2 */

+etanh( x2, x2 );	/* tanh( x/2 )			*/

+eadd( x2, eone, num );	/* 1 + tanh			*/

+eneg( x2 );

+eadd( x2, eone, den );	/* 1 - tanh			*/

+ediv( den, num, y );	/* (1 + tanh)/(1 - tanh)	*/

+

+/*y[NE-1] += i;*/

+if( sign )

+	{

+	ediv( y, eone, y );

+	i = -i;

+	}

+eldexp( y, i, y );	/* multiply by 2**i */

+}

diff --git a/test/math/ehead.h b/test/math/ehead.h
new file mode 100644
index 000000000..746b60df7
--- /dev/null
+++ b/test/math/ehead.h
@@ -0,0 +1,42 @@
+

+/* Include file for extended precision arithmetic programs.

+ */

+

+/* Number of 16 bit words in external x type format */

+#define NE 6

+

+/* Number of 16 bit words in internal format */

+#define NI (NE+3)

+

+/* Array offset to exponent */

+#define E 1

+

+/* Array offset to high guard word */

+#define M 2

+

+/* Number of bits of precision */

+#define NBITS ((NI-4)*16)

+

+/* Maximum number of decimal digits in ASCII conversion

+ * = NBITS*log10(2)

+ */

+#define NDEC (NBITS*8/27)

+

+/* The exponent of 1.0 */

+#define EXONE (0x3fff)

+

+void eadd(), esub(), emul(), ediv();

+int ecmp(), enormlz(), eshift();

+void eshup1(), eshup8(), eshup6(), eshdn1(), eshdn8(), eshdn6();

+void eabs(), eneg(), emov(), eclear(), einfin(), efloor();

+void eldexp(), efrexp(), eifrac(), ltoe();

+void esqrt(), elog(), eexp(), etanh(), epow();

+void asctoe(), asctoe24(), asctoe53(), asctoe64();

+void etoasc(), e24toasc(), e53toasc(), e64toasc();

+void etoe64(), etoe53(), etoe24(), e64toe(), e53toe(), e24toe();

+void mtherr();

+extern unsigned short ezero[], ehalf[], eone[], etwo[];

+extern unsigned short elog2[], esqrt2[];

+

+

+/* by Stephen L. Moshier. */

diff --git a/test/math/elog.c b/test/math/elog.c
new file mode 100644
index 000000000..8afc1e5b4
--- /dev/null
+++ b/test/math/elog.c
@@ -0,0 +1,92 @@
+/*						xlog.c	*/

+/* natural logarithm */

+/* by Stephen L. Moshier. */

+

+#include "mconf.h"

+#include "ehead.h"

+

+

+

+void elog( x, y )

+unsigned short *x, *y;

+{

+unsigned short xx[NE], z[NE], a[NE], b[NE], t[NE], qj[NE];

+long ex;

+int fex;

+

+

+if( x[NE-1] & (unsigned short )0x8000 )

+	{

+	eclear(y);

+	mtherr( "elog", DOMAIN );

+	return;

+	}

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

+	{

+	einfin( y );

+	eneg(y);

+	mtherr( "elog", SING );

+	return;

+	}

+if( ecmp( x, eone ) == 0 )

+	{

+	eclear( y );

+	return;

+	}

+

+/* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */

+efrexp( x, &fex, xx );

+/*

+emov(x, xx );

+ex = xx[NX-1] & 0x7fff;

+ex -= 0x3ffe;

+xx[NX-1] = 0x3ffe;

+*/

+

+/* Adjust range to 1/sqrt(2), sqrt(2) */

+esqrt2[NE-1] -= 1;

+if( ecmp( xx, esqrt2 ) < 0 )

+	{

+	fex -= 1;

+	emul( xx, etwo, xx );

+	}

+esqrt2[NE-1] += 1;

+

+esub( eone, xx, a );

+if( a[NE-1] == 0 )

+	{

+	eclear( y );

+	goto logdon;

+	}

+eadd( eone, xx, b );

+ediv( b, a, y );	/* store (x-1)/(x+1) in y */

+

+emul( y, y, z );

+

+emov( eone, a );

+emov( eone, b );

+emov( eone, qj );

+do

+	{

+	eadd( etwo, qj, qj );	/* 2 * i + 1		*/

+	emul( z, a, a );

+	ediv( qj, a, t );

+	eadd( t, b, b );

+	}

+while( ((b[NE-1] & 0x7fff) - (t[NE-1] & 0x7fff)) < NBITS );

+

+

+emul( b, y, y );

+emul( y, etwo, y );

+

+logdon:

+

+/* now add log of 2**ex */

+if( fex != 0 )

+	{

+	ex = fex;

+	ltoe( &ex, b );

+	emul( elog2, b, b );

+	eadd( b, y, y );

+	}

+}

diff --git a/test/math/eparanoi.c b/test/math/eparanoi.c
new file mode 100644
index 000000000..0e479f9f5
--- /dev/null
+++ b/test/math/eparanoi.c
@@ -0,0 +1,3550 @@
+/* paranoia.c arithmetic tester

+ *

+ * This is an implementation of the PARANOIA program.  It substitutes

+ * subroutine calls for ALL floating point arithmetic operations.

+ * This permits you to substitute your own experimental versions of

+ * arithmetic routines.  It also defeats compiler optimizations,

+ * so for native arithmetic you can be pretty sure you are testing

+ * the arithmetic and not the compiler.

+ *

+ * This version of PARANOIA omits the display of division by zero.

+ * It also omits the test for extra precise subexpressions, since

+ * they cannot occur in this context.  Otherwise it includes all the

+ * tests of the 27 Jan 86 distribution, plus a few additional tests.

+ * Commentary has been reduced to a minimum in order to make the program

+ * smaller.

+ *

+ * The original PARANOIA program, written by W. Kahan, C version

+ * by Thos Sumner and David Gay, can be downloaded free from the

+ * Internet NETLIB.  An MSDOS disk can be obtained for $15 from:

+ *   Richard Karpinski

+ *   6521 Raymond Street

+ *   Oakland, CA 94609

+ *

+ * Steve Moshier, 28 Oct 88

+ * last rev: 23 May 92

+ */

+

+#define DEBUG 0

+

+/* To use the native arithmetic of the computer, define NATIVE

+ * to be 1.  To use your own supplied arithmetic routines, NATIVE is 0.

+ */

+#define NATIVE 0

+

+/* gcc real.c interface */

+#define L128DOUBLE 0

+

+#include <stdio.h>

+

+

+

+

+/* Data structure of floating point number.  If NATIVE was

+ * selected above, you can define LDOUBLE 1 to test 80-bit long double

+ * precision or define it 0 to test 64-bit double precision.

+*/

+#define LDOUBLE 0

+#if NATIVE

+

+#define NE 1

+#if LDOUBLE

+#define FSIZE long double

+#define FLOAT(x) FSIZE x[NE]

+static FSIZE eone[NE] = {1.0L};	/* The constant 1.0 */

+#define ZSQRT sqrtl

+#define ZLOG logl

+#define ZFLOOR floorl

+#define ZPOW powl

+long double sqrtl(), logl(), floorl(), powl();

+#define FSETUP einit

+#else /* not LDOUBLE */

+#define FSIZE double

+#define FLOAT(x) FSIZE x[NE]

+static FSIZE eone[NE] = {1.0};	/* The constant 1.0 */

+#define ZSQRT sqrt

+#define ZLOG log

+#define ZFLOOR floor

+#define ZPOW pow

+double sqrt(), log(), floor(), pow();

+/* Coprocessor initialization,

+ * defeat underflow trap or what have you.

+ * This is required mainly on i386 and 68K processors.

+ */

+#define FSETUP dprec

+#endif /* double, not LDOUBLE */

+

+#else /* not NATIVE */

+

+/* Setup for extended double type.

+ * Put NE = 10 for real.c operating with TFmode support (16-byte reals)

+ * Put NE = 6 for real.c operating with XFmode support (10- or 12-byte reals)

+ * The value of NE must agree with that in ehead.h, if ieee.c is used.

+ */

+#define NE 6

+#define FSIZE unsigned short

+#define FLOAT(x) unsigned short x[NE]

+extern unsigned short eone[];

+#define FSETUP einit

+

+/* default for FSETUP */

+/*

+einit()

+{}

+*/

+

+error(s)

+char *s;

+{

+printf( "error: %s\n", s );

+}

+

+#endif	/* not NATIVE */

+

+

+

+#if L128DOUBLE

+/* real.c interface */

+

+#undef FSETUP

+#define FSETUP efsetup

+

+FLOAT(enone);

+

+#define ONE enone

+

+/* Use emov to convert from widest type to widest type, ... */

+/*

+#define ENTOE emov

+#define ETOEN emov

+*/

+

+/*                 ... else choose e24toe, e53toe, etc. */

+#define ENTOE e64toe

+#define ETOEN etoe64

+#define NNBITS 64

+

+#define NIBITS ((NE-1)*16)

+extern int rndprc;

+

+efsetup()

+{

+rndprc = NNBITS;

+ETOEN(eone, enone);

+}

+

+add(a,b,c)

+FLOAT(a);

+FLOAT(b);

+FLOAT(c);

+{

+unsigned short aa[10], bb[10], cc[10];

+

+ENTOE(a,aa);

+ENTOE(b,bb);

+eadd(aa,bb,cc);

+ETOEN(cc,c);

+}

+

+sub(a,b,c)

+FLOAT(a);

+FLOAT(b);

+FLOAT(c);

+{

+unsigned short aa[10], bb[10], cc[10];

+

+ENTOE(a,aa);

+ENTOE(b,bb);

+esub(aa,bb,cc);

+ETOEN(cc,c);

+}

+

+mul(a,b,c)

+FLOAT(a);

+FLOAT(b);

+FLOAT(c);

+{

+unsigned short aa[10], bb[10], cc[10];

+

+ENTOE(a,aa);

+ENTOE(b,bb);

+emul(aa,bb,cc);

+ETOEN(cc,c);

+}

+

+div(a,b,c)

+FLOAT(a);

+FLOAT(b);

+FLOAT(c);

+{

+unsigned short aa[10], bb[10], cc[10];

+

+ENTOE(a,aa);

+ENTOE(b,bb);

+ediv(aa,bb,cc);

+ETOEN(cc,c);

+}

+

+int cmp(a,b)

+FLOAT(a);

+FLOAT(b);

+{

+unsigned short aa[10], bb[10];

+int c;

+int ecmp();

+

+ENTOE(a,aa);

+ENTOE(b,bb);

+c = ecmp(aa,bb);

+return(c);

+}

+

+mov(a,b)

+FLOAT(a);

+FLOAT(b);

+{

+int i;

+

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

+	b[i] = a[i];

+}

+

+

+neg(a)

+FLOAT(a);

+{

+unsigned short aa[10];

+

+ENTOE(a,aa);

+eneg(aa);

+ETOEN(aa,a);

+}

+

+clear(a)

+FLOAT(a);

+{

+int i;

+

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

+	a[i] = 0;

+}

+

+FABS(a)

+FLOAT(a);

+{

+unsigned short aa[10];

+

+ENTOE(a,aa);

+eabs(aa);

+ETOEN(aa,a);

+}

+

+FLOOR(a,b)

+FLOAT(a);

+FLOAT(b);

+{

+unsigned short aa[10], bb[10];

+

+ENTOE(a,aa);

+efloor(aa,bb);

+ETOEN(bb,b);

+}

+

+LOG(a,b)

+FLOAT(a);

+FLOAT(b);

+{

+unsigned short aa[10], bb[10];

+int rndsav;

+

+ENTOE(a,aa);

+rndsav = rndprc;

+rndprc = NIBITS;

+elog(aa,bb);

+rndprc = rndsav;

+ETOEN(bb,b);

+}

+

+POW(a,b,c)

+FLOAT(a);

+FLOAT(b);

+FLOAT(c);

+{

+unsigned short aa[10], bb[10], cc[10];

+int rndsav;

+

+ENTOE(a,aa);

+ENTOE(b,bb);

+rndsav = rndprc;

+rndprc = NIBITS;

+epow(aa,bb,cc);

+rndprc = rndsav;

+ETOEN(cc,c);

+}

+

+SQRT(a,b)

+FLOAT(a);

+FLOAT(b);

+{

+unsigned short aa[10], bb[10];

+

+ENTOE(a,aa);

+esqrt(aa,bb);

+ETOEN(bb,b);

+}

+

+FTOL(x,ip,f)

+FLOAT(x);

+long *ip;

+FLOAT(f);

+{

+unsigned short xx[10], ff[10];

+

+ENTOE(x,xx);

+eifrac(xx,ip,ff);

+ETOEN(ff,f);

+}

+

+LTOF(ip,x)

+long *ip;

+FLOAT(x);

+{

+unsigned short xx[10];

+ltoe(ip,xx);

+ETOEN(xx,x);

+}

+

+TOASC(a,b,c)

+FLOAT(a);

+int b;

+char *c;

+{

+unsigned short xx[10];

+

+ENTOE(a,xx);

+etoasc(xx,b,c);

+}

+

+#else /* not L128DOUBLE */

+

+#define ONE eone

+

+/* Note all arguments of operation subroutines are pointers. */

+/* c = b + a */

+#define add(a,b,c) eadd(a,b,c)

+/* c = b - a */

+#define sub(a,b,c) esub(a,b,c)

+/* c = b * a */

+#define mul(a,b,c) emul(a,b,c)

+/* c = b / a */

+#define div(a,b,c) ediv(a,b,c)

+/* 1 if a>b, 0 if a==b, -1 if a<b */

+#define cmp(a,b) ecmp(a,b)

+/* b = a */

+#define mov(a,b) emov(a,b)

+/* a = -a */

+#define neg(a) eneg(a)

+/* a = 0 */

+#define clear(a) eclear(a)

+

+#define FABS(x) eabs(x)

+#define FLOOR(x,y) efloor(x,y)

+#define LOG(x,y) elog(x,y)

+#define POW(x,y,z) epow(x,y,z)

+#define SQRT(x,y) esqrt(x,y)

+

+/* x = &FLOAT input, i = &long integer part, f = &FLOAT fractional part */

+#define FTOL(x,i,f) eifrac(x,i,f)

+

+/* i = &long integer input, x = &FLOAT output */

+#define LTOF(i,x) ltoe(i,x)

+

+/* Convert FLOAT a to decimal ASCII string with b digits */

+#define TOASC(a,b,c) etoasc(a,b,c)

+#endif /* not L128DOUBLE */

+

+

+

+/* The following subroutines are implementations of the above

+ * named functions, using the native or default arithmetic.

+ */

+#if NATIVE

+eadd(a,b,c)

+FSIZE *a, *b, *c;

+{

+*c = *b + *a;

+}

+

+esub(a,b,c)

+FSIZE *a, *b, *c;

+{

+*c = *b - *a;

+}

+

+emul(a,b,c)

+FSIZE *a, *b, *c;

+{

+*c = (*b) * (*a);

+}

+

+ediv(a,b,c)

+FSIZE *a, *b, *c;

+{

+*c = (*b) / (*a);

+}

+

+

+/* Important note: comparison can be done by subracting

+ * or by a compare instruction that may or may not be

+ * equivalent to subtracting.

+ */

+ecmp(a,b)

+FSIZE *a, *b;

+{

+if( (*a) > (*b) )

+	return( 1 );

+if( (*a) < (*b) )

+	return( -1 );

+if( (*a) != (*b) )

+	goto cmpf;

+if( (*a) == (*b) )

+	return( 0 );

+cmpf:

+printf( "Compare fails\n" );

+return(0);

+}

+

+

+emov( a, b )

+FSIZE *a, *b;

+{

+*b = *a;

+}

+

+eneg( a )

+FSIZE *a;

+{

+*a = -(*a);

+}

+

+eclear(a)

+FSIZE *a;

+{

+*a = 0.0;

+}

+

+eabs(x)

+FSIZE *x;

+{

+if( (*x) < 0.0 )

+	*x = -(*x);

+}

+

+efloor(x,y)

+FSIZE *x, *y;

+{

+

+*y = (FSIZE )ZFLOOR( *x );

+}

+

+elog(x,y)

+FSIZE *x, *y;

+{

+

+*y = (FSIZE )ZLOG( *x );

+}

+

+epow(x,y,z)

+FSIZE *x, *y, *z;

+{

+

+*z = (FSIZE )ZPOW( *x, *y );

+}

+

+esqrt(x,y)

+FSIZE *x, *y;

+{

+

+*y = (FSIZE )ZSQRT( *x );

+}

+

+

+eifrac(x,i,f)

+FSIZE *x;

+long *i;

+FSIZE *f;

+{

+FSIZE y;

+

+y = (FSIZE )ZFLOOR( *x );

+if( y < 0.0 )

+	{

+	*f = y - *x;

+	*i = -y;

+	}

+else

+	{

+	*f = *x - y;

+	*i = y;

+	}

+}

+

+

+ltoe(i,x)

+long *i;

+FSIZE *x;

+{

+*x = *i;

+}

+

+

+etoasc(a,str,n)

+FSIZE *a;

+char *str;

+int n;

+{

+double x;

+

+x = (double )(*a);

+sprintf( str, " %.17e ", x );

+}

+

+/* default for FSETUP */

+einit()

+{}

+

+#endif	/* NATIVE */

+

+

+

+

+FLOAT(Radix);

+FLOAT(BInvrse);

+FLOAT(RadixD2);

+FLOAT(BMinusU2);

+/*Small floating point constants.*/

+FLOAT(Zero);

+FLOAT(Half);

+FLOAT(One);

+FLOAT(Two);

+FLOAT(Three);

+FLOAT(Four);

+FLOAT(Five);

+FLOAT(Six);

+FLOAT(Eight);

+FLOAT(Nine);

+FLOAT(Ten);

+FLOAT(TwentySeven);

+FLOAT(ThirtyTwo);

+FLOAT(TwoForty);

+FLOAT(MinusOne );

+FLOAT(OneAndHalf);

+

+/*Integer constants*/

+int NoTrials = 20; /*Number of tests for commutativity. */

+#define False 0

+#define True 1

+

+/* Definitions for declared types 

+	Guard == (Yes, No);

+	Rounding == (Chopped, Rounded, Other);

+	Message == packed array [1..40] of char;

+	Class == (Flaw, Defect, Serious, Failure);

+	  */

+#define Yes 1

+#define No  0

+#define Chopped 2

+#define Rounded 1

+#define Other   0

+#define Flaw    3

+#define Defect  2

+#define Serious 1

+#define Failure 0

+

+typedef int Guard, Rounding, Class;

+typedef char Message;

+

+/* Declarations of Variables */

+FLOAT(AInvrse);

+FLOAT(A1);

+FLOAT(C);

+FLOAT(CInvrse);

+FLOAT(D);

+FLOAT(FourD);

+FLOAT(E0);

+FLOAT(E1);

+FLOAT(Exp2);

+FLOAT(E3);

+FLOAT(MinSqEr);

+FLOAT(SqEr);

+FLOAT(MaxSqEr);

+FLOAT(E9);

+FLOAT(Third);

+FLOAT(F6);

+FLOAT(F9);

+FLOAT(H);

+FLOAT(HInvrse);

+FLOAT(StickyBit);

+FLOAT(J);

+FLOAT(MyZero);

+FLOAT(Precision);

+FLOAT(Q);

+FLOAT(Q9);

+FLOAT(R);

+FLOAT(Random9);

+FLOAT(T);

+FLOAT(Underflow);

+FLOAT(S);

+FLOAT(OneUlp);

+FLOAT(UfThold);

+FLOAT(U1);

+FLOAT(U2);

+FLOAT(V);

+FLOAT(V0);

+FLOAT(V9);

+FLOAT(W);

+FLOAT(X);

+FLOAT(X1);

+FLOAT(X2);

+FLOAT(X8);

+FLOAT(Random1);

+FLOAT(Y);

+FLOAT(YY1);

+FLOAT(Y2);

+FLOAT(Random2);

+FLOAT(Z);

+FLOAT(PseudoZero);

+FLOAT(Z1);

+FLOAT(Z2);

+FLOAT(Z9);

+static FLOAT(t);

+FLOAT(t2);

+FLOAT(Sqarg);

+int ErrCnt[4];

+int fpecount;

+int Milestone;

+int PageNo;

+int I, M, N, N1, stkflg;

+Guard GMult, GDiv, GAddSub;

+Rounding RMult, RDiv, RAddSub, RSqrt;

+int Break, Done, NotMonot, Monot, Anomaly, IEEE;

+int SqRWrng, UfNGrad;

+int k, k2;

+int Indx;

+char ch[8];

+

+long lngint, lng2; /* intermediate for conversion between int and FLOAT */

+

+/* Computed constants. */

+/*U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */

+/*U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */

+

+

+show( x )

+short x[];

+{

+int i;

+char s[80];

+

+/* Number of 16-bit groups to display */

+#if NATIVE

+#if LDOUBLE

+#define NPRT (sizeof( long double )/2)

+#else

+#define NPRT (sizeof( double )/2)

+#endif

+#else

+#define NPRT NE

+#endif

+

+TOASC( x, s, 70 );

+printf( "%s\n", s );

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

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

+printf( "\n" );

+}

+

+/* define NOSIGNAL */

+#ifndef NOSIGNAL

+#include <signal.h>

+#endif

+#include <setjmp.h>

+jmp_buf ovfl_buf;

+/*typedef int (*Sig_type)();*/

+typedef void (*Sig_type)();

+Sig_type sigsave;

+

+/* Fl