Rework the math library tests per the glibc math test code, with

many unsupported tests disabled for the moment.
 -Erik

24 files changed, 5522 insertions, 9833 deletions

diff --git a/test/math/.cvsignore b/test/math/.cvsignore
new file mode 100644
index 000000000..31ed5f055
--- /dev/null
+++ b/test/math/.cvsignore
@@ -0,0 +1,8 @@
+libm-test-ulps.h
+libm-test.c
+test-double
+test-idouble
+test-float
+test-ifloat
+test-ldouble
+test-ildouble
diff --git a/test/math/Makefile b/test/math/Makefile
index 4d0bc5ee6..d41074939 100644
--- a/test/math/Makefile
+++ b/test/math/Makefile
@@ -18,74 +18,57 @@
 
 
 
-# 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
+CFLAGS+=-D_GNU_SOURCE -DNO_LONG_DOUBLE
+EXTRA_LIBS=-lm
+PERL=/usr/bin/perl
+
+TARGETS:=
+libm-tests:=
+libm-tests+= test-double test-idouble
+#libm-tests+= test-float test-ifloat
+#libm-tests+= test-ldouble test-ildouble
+libm-tests.o = $(addsuffix .o,$(libm-tests))
+
+libm-tests-generated = libm-test-ulps.h libm-test.c
+generated += $(libm-tests-generated) libm-test.stmp
+TARGETS += $(libm-tests) #$(libm-tests-generated)
+
+all: libm-test.c $(TARGETS)
+	
+test-double: test-double.o
+	$(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS)
+	-./$@
+test-idouble: test-idouble.o
+	$(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS)
+	-./$@
+test-float: test-float.o
+	$(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS)
+	-./$@
+test-ifloat: test-ifloat.o
+	$(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS)
+	-./$@
+test-ldouble: test-ldouble.o
+	$(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS)
+	-./$@
+test-ildouble: test-ildoubl.o
+	$(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS)
+	-./$@
+
+test-float.o: libm-test.c
+test-ifloat.o: libm-test.c 
+test-double.o: libm-test.c
+test-idouble.o: libm-test.c
+test-ldouble.o: libm-test.c
+test-ildoubl.o: libm-test.c
+
+ulps-file = $(firstword $(wildcard $(config-sysdirs:%=$(..)%/libm-test-ulps)))
+
+libm-test.c: $(ulps-file) libm-test.inc gen-libm-test.pl
+	$(PERL) ./gen-libm-test.pl -u $< ./libm-test.inc -o "." 2>&1 > /dev/null
 
 clean:
-	rm -f *.[oa] *~ core $(TARGETS)
+	rm -f *.[oa] *~ core $(TARGETS) $(generated)
 
 
diff --git a/test/math/drand.c b/test/math/drand.c
deleted file mode 100644
index 9eedf71fc..000000000
--- a/test/math/drand.c
+++ /dev/null
@@ -1,158 +0,0 @@
-/*							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
deleted file mode 100644
index cfddbe3e2..000000000
--- a/test/math/econst.c
+++ /dev/null
@@ -1,96 +0,0 @@
-/*							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
deleted file mode 100644
index 14ea9899d..000000000
--- a/test/math/eexp.c
+++ /dev/null
@@ -1,77 +0,0 @@
-/*							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
deleted file mode 100644
index 24c95ce05..000000000
--- a/test/math/ehead.h
+++ /dev/null
@@ -1,42 +0,0 @@
-
-/* 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
deleted file mode 100644
index bc517b197..000000000
--- a/test/math/elog.c
+++ /dev/null
@@ -1,92 +0,0 @@
-/*						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
deleted file mode 100644
index 84cab73f8..000000000
--- a/test/math/eparanoi.c
+++ /dev/null
@@ -1,3550 +0,0 @@
-/* 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);