From f108799afa9b1d207e7139608c261f5546eb56bf Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Tue, 2 Oct 2001 10:45:16 +0000 Subject: Add in some math lib tests --- test/math/Makefile | 71 + test/math/drand.c | 158 ++ test/math/econst.c | 96 ++ test/math/eexp.c | 77 + test/math/ehead.h | 42 + test/math/elog.c | 92 ++ test/math/eparanoi.c | 3550 +++++++++++++++++++++++++++++++++++++++++++ test/math/epow.c | 215 +++ test/math/etanh.c | 52 + test/math/etodec.c | 181 +++ test/math/ieee.c | 4119 ++++++++++++++++++++++++++++++++++++++++++++++++++ test/math/ieetst.c | 850 +++++++++++ test/math/ieetst.doc | 132 ++ test/math/mconf.h | 108 ++ test/math/mtherr.c | 96 ++ 15 files changed, 9839 insertions(+) create mode 100644 test/math/Makefile create mode 100644 test/math/drand.c create mode 100644 test/math/econst.c create mode 100644 test/math/eexp.c create mode 100644 test/math/ehead.h create mode 100644 test/math/elog.c create mode 100644 test/math/eparanoi.c create mode 100644 test/math/epow.c create mode 100644 test/math/etanh.c create mode 100644 test/math/etodec.c create mode 100644 test/math/ieee.c create mode 100644 test/math/ieetst.c create mode 100644 test/math/ieetst.doc create mode 100644 test/math/mconf.h create mode 100644 test/math/mtherr.c (limited to 'test') 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 + + + + +/* 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; ib, 0 if a==b, -1 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 +#endif +#include +jmp_buf ovfl_buf; +/*typedef int (*Sig_type)();*/ +typedef void (*Sig_type)(); +Sig_type sigsave; + +/* Floating point exception receiver */ +void sigfpe() +{ +fpecount++; +printf( "\n* * * FLOATING-POINT ERROR * * *\n" ); +/* reinitialize the floating point unit */ +FSETUP(); +fflush(stdout); +if( sigsave ) + { +#ifndef NOSIGNAL + signal( SIGFPE, sigsave ); +#endif + sigsave = 0; + longjmp( ovfl_buf, 1 ); + } +abort(); +} + + +main() +{ + +/* Do coprocessor or other initializations */ +FSETUP(); + +printf( + "This version of paranoia omits test for extra precise subexpressions\n" ); +printf( "and includes a few additional tests.\n" ); + +clear(Zero); +printf( "0 = " ); +show( Zero ); +mov( ONE, One); +printf( "1 = " ); +show( One ); +add( One, One, Two ); +printf( "1+1 = " ); +show( Two ); +add( Two, One, Three ); +add( Three, One, Four ); +add( Four, One, Five ); +add( Five, One, Six ); +add( Four, Four, Eight ); +mul( Three, Three, Nine ); +add( Nine, One, Ten ); +mul( Nine, Three, TwentySeven ); +mul( Four, Eight, ThirtyTwo ); +mul( Four, Five, t ); +mul( t, Three, t ); +mul( t, Four, TwoForty ); +mov( One, MinusOne ); +neg( MinusOne ); +div( Two, One, Half ); +add( One, Half, OneAndHalf ); +ErrCnt[Failure] = 0; +ErrCnt[Serious] = 0; +ErrCnt[Defect] = 0; +ErrCnt[Flaw] = 0; +PageNo = 1; +#ifndef NOSIGNAL +signal( SIGFPE, sigfpe ); +#endif +printf("Program is now RUNNING tests on small integers:\n"); + +add( Zero, Zero, t ); +if( cmp( t, Zero ) != 0) + { + printf( "0+0 != 0\n" ); + ErrCnt[Failure] += 1; + } +sub( One, One, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "1-1 != 0\n" ); + ErrCnt[Failure] += 1; + } +if( cmp( One, Zero ) <= 0 ) + { + printf( "1 <= 0\n" ); + ErrCnt[Failure] += 1; + } +add( One, One, t ); +if( cmp( t, Two ) != 0 ) + { + printf( "1+1 != 2\n" ); + ErrCnt[Failure] += 1; + } +mov( Zero, Z ); +neg( Z ); +FLOOR( Z, t ); +if( cmp(t,Zero) != 0 ) + { + ErrCnt[Serious] += 1; + printf( "FLOOR(-0) should equal 0, is = " ); + show( t ); + } +if( cmp(Z, Zero) != 0) + { + ErrCnt[Failure] += 1; + printf("Comparison alleges that -0.0 is Non-zero!\n"); + } +else + { + div( TwoForty, One, U1 ); /* U1 = 0.001 */ + mov( One, Radix ); + TstPtUf(); + } +add( Two, One, t ); +if( cmp( t, Three ) != 0 ) + { + printf( "2+1 != 3\n" ); + ErrCnt[Failure] += 1; + } +add( Three, One, t ); +if( cmp( t, Four ) != 0 ) + { + printf( "3+1 != 4\n" ); + ErrCnt[Failure] += 1; + } +mov( Two, t ); +neg( t ); +mul( Two, t, t ); +add( Four, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "4+2*(-2) != 0\n" ); + ErrCnt[Failure] += 1; + } +sub( Three, Four, t ); +sub( One, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "4-3-1 != 0\n" ); + ErrCnt[Failure] += 1; + } + sub( One, Zero, t ); +if( cmp( t, MinusOne ) != 0 ) + { + printf( "-1 != 0-1\n" ); + ErrCnt[Failure] += 1; + } +add( One, MinusOne, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "1+(-1) != 0\n" ); + ErrCnt[Failure] += 1; + } +mov( One, t ); +FABS( t ); +add( MinusOne, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "-1+abs(1) != 0\n" ); + ErrCnt[Failure] += 1; + } +mul( MinusOne, MinusOne, t ); +add( MinusOne, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "-1+(-1)*(-1) != 0\n" ); + ErrCnt[Failure] += 1; + } +add( Half, MinusOne, t ); +add( Half, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "1/2 + (-1) + 1/2 != 0\n" ); + ErrCnt[Failure] += 1; + } +Milestone = 10; +mul( Three, Three, t ); +if( cmp( t, Nine ) != 0 ) + { + printf( "3*3 != 9\n" ); + ErrCnt[Failure] += 1; + } +mul( Nine, Three, t ); +if( cmp( t, TwentySeven ) != 0 ) + { + printf( "3*9 != 27\n" ); + ErrCnt[Failure] += 1; + } +add( Four, Four, t ); +if( cmp( t, Eight ) != 0 ) + { + printf( "4+4 != 8\n" ); + ErrCnt[Failure] += 1; + } +mul( Eight, Four, t ); +if( cmp( t, ThirtyTwo ) != 0 ) + { + printf( "8*4 != 32\n" ); + ErrCnt[Failure] += 1; + } +sub( TwentySeven, ThirtyTwo, t ); +sub( Four, t, t ); +sub( One, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "32-27-4-1 != 0\n" ); + ErrCnt[Failure] += 1; + } +add( Four, One, t ); +if( cmp( t, Five ) != 0 ) + { + printf( "4+1 != 5\n" ); + ErrCnt[Failure] += 1; + } +mul( Four, Five, t ); +mul( Three, t, t ); +mul( Four, t, t ); +if( cmp( t, TwoForty ) != 0 ) + { + printf( "4*5*3*4 != 240\n" ); + ErrCnt[Failure] += 1; + } +div( Three, TwoForty, t ); +mul( Four, Four, t2 ); +mul( Five, t2, t2 ); +sub( t2, t2, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "240/3 - 4*4*5 != 0\n" ); + ErrCnt[Failure] += 1; + } +div( Four, TwoForty, t ); +mul( Five, Three, t2 ); +mul( Four, t2, t2 ); +sub( t2, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "240/4 - 5*3*4 != 0\n" ); + ErrCnt[Failure] += 1; + } +div( Five, TwoForty, t ); +mul( Four, Three, t2 ); +mul( Four, t2, t2 ); +sub( t2, t, t ); +if( cmp( t, Zero ) != 0 ) + { + printf( "240/5 - 4*3*4 != 0\n" ); + ErrCnt[Failure] += 1; + } +if(ErrCnt[Failure] == 0) + { +printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n\n"); + } +printf("Searching for Radix and Precision.\n"); +mov( One, W ); +do + { + add( W, W, W ); + add( W, One, Y ); + sub( W, Y, Z ); + sub( One, Z, Y ); + mov( Y, t ); + FABS(t); + add( MinusOne, t, t ); + k = cmp( t, Zero ); + } +while( k < 0 ); +/*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/ +mov( Zero, Precision ); +mov( One, Y ); +do + { + add( W, Y, Radix ); + add( Y, Y, Y ); + sub( W, Radix, Radix ); + k = cmp( Radix, Zero ); + } +while( k == 0); + +if( cmp(Radix, Two) < 0 ) + mov( One, Radix ); +printf("Radix = " ); +show( Radix ); +if( cmp(Radix, One) != 0) + { + mov( One, W ); + do + { + add( One, Precision, Precision ); + mul( W, Radix, W ); + add( W, One, Y ); + sub( W, Y, t ); + k = cmp( t, One ); + } + while( k == 0 ); + } +/* now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 */ +div( W, One, U1 ); +mul( Radix, U1, U2 ); +printf( "Closest relative separation found is U 1 = " ); +show( U1 ); +printf( "Recalculating radix and precision." ); + +/*save old values*/ +mov( Radix, E0 ); +mov( U1, E1 ); +mov( U2, E9 ); +mov( Precision, E3 ); + +div( Three, Four, X ); +sub( One, X, Third ); +sub( Third, Half, F6 ); +add( F6, F6, X ); +sub( Third, X, X ); +FABS( X ); +if( cmp(X, U2) < 0 ) + mov( U2, X ); + +/*... now X = (unknown no.) ulps of 1+...*/ +do + { + mov( X, U2 ); +/* Y = Half * U2 + ThirtyTwo * U2 * U2; */ + mul( ThirtyTwo, U2, t ); + mul( t, U2, t ); + mul( Half, U2, Y ); + add( t, Y, Y ); + add( One, Y, Y ); + sub( One, Y, X ); + k = cmp( U2, X ); + k2 = cmp( X, Zero ); + } +while ( ! ((k <= 0) || (k2 <= 0))); + +/*... now U2 == 1 ulp of 1 + ... */ +div( Three, Two, X ); +sub( Half, X, F6 ); +add( F6, F6, Third ); +sub( Half, Third, X ); +add( F6, X, X ); +FABS( X ); +if( cmp(X, U1) < 0 ) + mov( U1, X ); + +/*... now X == (unknown no.) ulps of 1 -... */ +do + { + mov( X, U1 ); + /* Y = Half * U1 + ThirtyTwo * U1 * U1;*/ + mul( ThirtyTwo, U1, t ); + mul( U1, t, t ); + mul( Half, U1, Y ); + add( t, Y, Y ); + sub( Y, Half, Y ); + add( Half, Y, X ); + sub( X, Half, Y ); + add( Half, Y, X ); + k = cmp( U1, X ); + k2 = cmp( X, Zero ); + } while ( ! ((k <= 0) || (k2 <= 0))); +/*... now U1 == 1 ulp of 1 - ... */ +if( cmp( U1, E1 ) == 0 ) + printf("confirms closest relative separation U1 .\n"); +else + { + printf("gets better closest relative separation U1 = " ); + show( U1 ); + } +div( U1, One, W ); +sub( U1, Half, F9 ); +add( F9, Half, F9 ); +div( U1, U2, t ); +div( TwoForty, One, t2 ); +add( t2, t, t ); +FLOOR( t, Radix ); +if( cmp(Radix, E0) == 0 ) + printf("Radix confirmed.\n"); +else + { + printf("MYSTERY: recalculated Radix = " ); + show( Radix ); + mov( E0, Radix ); + } +add( Eight, Eight, t ); +if( cmp( Radix, t ) > 0 ) + { + printf( "Radix is too big: roundoff problems\n" ); + ErrCnt[Defect] += 1; + } +k = 1; +if( cmp( Radix, Two ) == 0 ) + k = 0; +if( cmp( Radix, Ten ) == 0 ) + k = 0; +if( cmp( Radix, One ) == 0 ) + k = 0; +if( k != 0 ) + { + printf( "Radix is not as good as 2 or 10\n" ); + ErrCnt[Flaw] += 1; + } +/*=============================================*/ +Milestone = 20; +/*=============================================*/ +sub( Half, F9, t ); +if( cmp( t, Half ) >= 0 ) + { + printf( "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?\n" ); + ErrCnt[Failure] += 1; + } +mov( F9, X ); +I = 1; +sub( Half, X, Y ); +sub( Half, Y, Z ); +if( (cmp( X, One ) == 0) && (cmp( Z, Zero) != 0) ) + { + printf( "Comparison is fuzzy ,X=1 but X-1/2-1/2 != 0\n" ); + ErrCnt[Failure] += 1; + } +add( One, U2, X ); +I = 0; +/*=============================================*/ +Milestone = 25; +/*=============================================*/ +/*... BMinusU2 = nextafter(Radix, 0) */ + +sub( One, Radix, BMinusU2 ); +sub( U2, BMinusU2, t ); +add( One, t, BMinusU2 ); +/* Purify Integers */ +if( cmp(Radix,One) != 0 ) + { +/*X = - TwoForty * LOG(U1) / LOG(Radix);*/ + LOG( U1, X ); + LOG( Radix, t ); + div( t, X, X ); + mul( TwoForty, X, X ); + neg( X ); + + add( Half, X, Y ); + FLOOR( Y, Y ); + sub( Y, X, t ); + FABS( t ); + mul( Four, t, t ); + if( cmp( t, One ) < 0 ) + mov( Y, X ); + div( TwoForty, X, Precision ); + add( Half, Precision, Y ); + FLOOR( Y, Y ); + sub( Y, Precision, t ); + FABS( t ); + mul( TwoForty, t, t ); + if( cmp( t, Half ) < 0 ) + mov( Y, Precision ); + } +FLOOR( Precision, t ); +if( (cmp( Precision, t ) != 0) || (cmp( Radix, One ) == 0) ) + { + printf("Precision cannot be characterized by an Integer number\n"); + printf("of significant digits but, by itself, this is a minor flaw.\n"); + } +if( cmp(Radix, One) == 0 ) + printf("logarithmic encoding has precision characterized solely by U1.\n"); +else + { + printf("The number of significant digits of the Radix is " ); + show( Precision ); + } +mul( U2, Nine, t ); +mul( Nine, t, t ); +mul( TwoForty, t, t ); +if( cmp( t, One ) >= 0 ) + { + printf( "Precision worse than 5 decimal figures\n" ); + ErrCnt[Serious] += 1; + } +/*=============================================*/ +Milestone = 30; +/*=============================================*/ +/* Test for extra-precise subepressions has been deleted. */ +Milestone = 35; +/*=============================================*/ +if( cmp(Radix,Two) >= 0 ) + { + mul( Radix, Radix, t ); + div( t, W, X ); + add( X, One, Y ); + sub( X, Y, Z ); + add( Z, U2, T ); + sub( Z, T, X ); + if( cmp( X, U2 ) != 0 ) + { + printf( "Subtraction is not normalized X=Y,X+Z != Y+Z!\n" ); + ErrCnt[Failure] += 1; + } + if( cmp(X,U2) == 0 ) + printf("Subtraction appears to be normalized, as it should be."); + } + +printf("\nChecking for guard digit in *, /, and -.\n"); +mul( F9, One, Y ); +mul( One, F9, Z ); +sub( Half, F9, X ); +sub( Half, Y, Y ); +sub( X, Y, Y ); +sub( Half, Z, Z ); +sub( X, Z, Z ); +add( One, U2, X ); +mul( X, Radix, T ); +mul( Radix, X, R ); +sub( Radix, T, X ); +mul( Radix, U2, t ); +sub( t, X, X ); +sub( Radix, R, T ); +mul( Radix, U2, t ); +sub( t, T, T ); +sub( One, Radix, t ); +mul( t, X, X ); +sub( One, Radix, t ); +mul( t, T, T ); + +k = cmp(X,Zero); +k |= cmp(Y,Zero); +k |= cmp(Z,Zero); +k |= cmp(T,Zero); +if( k == 0 ) + GMult = Yes; +else + { + GMult = No; + ErrCnt[Serious] += 1; + printf( "* lacks a Guard Digit, so 1*X != X\n" ); + } +mul( Radix, U2, Z ); +add( One, Z, X ); +add( X, Z, Y ); +mul( X, X, t ); +sub( t, Y, Y ); +FABS( Y ); +sub( U2, Y, Y ); +sub( U2, One, X ); +sub( U2, X, Z ); +mul( X, X, t ); +sub( t, Z, Z ); +FABS( Z ); +sub( U1, Z, Z ); +if( (cmp(Y,Zero) > 0) || (cmp(Z,Zero) > 0) ) + { + ErrCnt[Failure] += 1; + printf( "* gets too many final digits wrong.\n" ); + } +sub( U2, One, Y ); +add( One, U2, X ); +div( Y, One, Z ); +sub( X, Z, Y ); +div( Three, One, X ); +div( Nine, Three, Z ); +sub( Z, X, X ); +div( TwentySeven, Nine, T ); +sub( T, Z, Z ); +k = cmp( X, Zero ); +k |= cmp( Y, Zero ); +k |= cmp( Z, Zero ); +if( k ) + { + ErrCnt[Defect] += 1; +printf( "Division lacks a Guard Digit, so error can exceed 1 ulp\n" ); +printf( "or 1/3 and 3/9 and 9/27 may disagree\n" ); + } +div( One, F9, Y ); +sub( Half, F9, X ); +sub( Half, Y, Y ); +sub( X, Y, Y ); +add( One, U2, X ); +div( One, X, T ); +sub( X, T, X ); +k = cmp( X, Zero ); +k |= cmp( Y, Zero ); +k |= cmp( Z, Zero ); +if( k == 0 ) + GDiv = Yes; +else + { + GDiv = No; + ErrCnt[Serious] += 1; + printf( "Division lacks a Guard Digit, so X/1 != X\n" ); + } +add( One, U2, X ); +div( X, One, X ); +sub( Half, X, Y ); +sub( Half, Y, Y ); +if( cmp(Y,Zero) >= 0 ) + { + ErrCnt[Serious] += 1; + printf( "Computed value of 1/1.000..1 >= 1\n" ); + } +sub( U2, One, X ); +mul( Radix, U2, Y ); +add( One, Y, Y ); +mul( X, Radix, Z ); +mul( Y, Radix, T ); +div( Radix, Z, R ); +div( Radix, T, StickyBit ); +sub( X, R, X ); +sub( Y, StickyBit, Y ); +k = cmp( X, Zero ); +k |= cmp( Y, Zero ); +if( k ) + { + ErrCnt[Failure] += 1; + printf( "* and/or / gets too many last digits wrong\n" ); + } +sub( U1, One, Y ); +sub( F9, One, X ); +sub( Y, One, Y ); +sub( U2, Radix, T ); +sub( BMinusU2, Radix, Z ); +sub( T, Radix, T ); +k = cmp( X, U1 ); +k |= cmp( Y, U1 ); +k |= cmp( Z, U2 ); +k |= cmp( T, U2 ); +if( k == 0 ) + GAddSub = Yes; +else + { + GAddSub = No; + ErrCnt[Serious] += 1; + printf( "- lacks Guard Digit, so cancellation is obscured\n" ); + } +sub( One, F9, t ); +if( (cmp(F9,One) != 0) && (cmp(t,Zero) >= 0) ) + { + ErrCnt[Serious] += 1; + printf("comparison alleges (1-U1) < 1 although\n"); + printf(" subtration yields (1-U1) - 1 = 0 , thereby vitiating\n"); + printf(" such precautions against division by zero as\n"); + printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n"); + } +if (GMult == Yes && GDiv == Yes && GAddSub == Yes) + printf(" *, /, and - appear to have guard digits, as they should.\n"); +/*=============================================*/ +Milestone = 40; +/*=============================================*/ +printf("Checking rounding on multiply, divide and add/subtract.\n"); +RMult = Other; +RDiv = Other; +RAddSub = Other; +div( Two, Radix, RadixD2 ); +mov( Two, A1 ); +Done = False; +do + { + mov( Radix, AInvrse ); + do + { + mov( AInvrse, X ); + div( A1, AInvrse, AInvrse ); + FLOOR( AInvrse, t ); + k = cmp( t, AInvrse ); + } + while( ! (k != 0 ) ); + k = cmp( X, One ); + k2 = cmp( A1, Three ); + Done = (k == 0) || (k2 > 0); + if(! Done) + add( Nine, One, A1 ); + } +while( ! (Done)); +if( cmp(X, One) == 0 ) + mov( Radix, A1 ); +div( A1, One, AInvrse ); +mov( A1, X ); +mov( AInvrse, Y ); +Done = False; +do + { + mul( X, Y, Z ); + sub( Half, Z, Z ); + if( cmp( Z, Half ) != 0 ) + { + ErrCnt[Failure] += 1; + printf( "X * (1/X) differs from 1\n" ); + } + k = cmp( X, Radix ); + Done = (k == 0); + mov( Radix, X ); + div( X, One, Y ); + } +while( ! (Done)); + +add( One, U2, Y2 ); +sub( U2, One, YY1 ); +sub( U2, OneAndHalf, X ); +add( OneAndHalf, U2, Y ); +sub( U2, X, Z ); +mul( Z, Y2, Z ); +mul( Y, YY1, T ); +sub( X, Z, Z ); +sub( X, T, T ); +mul( X, Y2, X ); +add( Y, U2, Y ); +mul( Y, YY1, Y ); +sub( OneAndHalf, X, X ); +sub( OneAndHalf, Y, Y ); +k = cmp( X, Zero ); +k |= cmp( Y, Zero ); +k |= cmp( Z, Zero ); +if( cmp( T, Zero ) > 0 ) + k = 1; +if( k == 0 ) + { + add( OneAndHalf, U2, X ); + mul( X, Y2, X ); + sub( U2, OneAndHalf, Y ); + sub( U2, Y, Y ); + add( OneAndHalf, U2, Z ); + add( U2, Z, Z ); + sub( U2, OneAndHalf, T ); + mul( T, YY1, T ); + add( Z, U2, t ); + sub( t, X, X ); + mul( Y, YY1, StickyBit ); + mul( Z, Y2, S ); + sub( Y, T, T ); + sub( Y, U2, Y ); + add( StickyBit, Y, Y ); +/* Z = S - (Z + U2 + U2); */ + add( Z, U2, t ); + add( t, U2, t ); + sub( t, S, Z ); + add( Y2, U2, t ); + mul( t, YY1, StickyBit ); + mul( Y2, YY1, YY1 ); + sub( Y2, StickyBit, StickyBit ); + sub( Half, YY1, YY1 ); + k = cmp( X, Zero ); + k |= cmp( Y, Zero ); + k |= cmp( Z, Zero ); + k |= cmp( T, Zero ); + k |= cmp( StickyBit, Zero ); + k |= cmp( YY1, Half ); + if( k == 0 ) + { + RMult = Rounded; + printf("Multiplication appears to round correctly.\n"); + } + else + { + add( X, U2, t ); + k = cmp( t, Zero ); + if( cmp( Y, Zero ) >= 0 ) + k |= 1; + add( Z, U2, t ); + k |= cmp( t, Zero ); + if( cmp( T, Zero ) >= 0 ) + k |= 1; + add( StickyBit, U2, t ); + k |= cmp( t, Zero ); + if( cmp(YY1, Half) >= 0 ) + k |= 1; + if( k == 0 ) + { + printf("Multiplication appears to chop.\n"); + } + else + { + printf("* is neither chopped nor correctly rounded.\n"); + } + if( (RMult == Rounded) && (GMult == No) ) + printf("Multiplication has inconsistent result"); + } + } +else + printf("* is neither chopped nor correctly rounded.\n"); + +/*=============================================*/ +Milestone = 45; +/*=============================================*/ +add( One, U2, Y2 ); +sub( U2, One, YY1 ); +add( OneAndHalf, U2, Z ); +add( Z, U2, Z ); +div( Y2, Z, X ); +sub( U2, OneAndHalf, T ); +sub( U2, T, T ); +sub( U2, T, Y ); +div( YY1, Y, Y ); +add( Z, U2, Z ); +div( Y2, Z, Z ); +sub( OneAndHalf, X, X ); +sub( T, Y, Y ); +div( YY1, T, T ); +add( OneAndHalf, U2, t ); +sub( t, Z, Z ); +sub( OneAndHalf, U2, t ); +add( t, T, T ); +k = 0; +if( cmp( X, Zero ) > 0 ) + k = 1; +if( cmp( Y, Zero ) > 0 ) + k = 1; +if( cmp( Z, Zero ) > 0 ) + k = 1; +if( cmp( T, Zero ) > 0 ) + k = 1; +if( k == 0 ) + { + div( Y2, OneAndHalf, X ); + sub( U2, OneAndHalf, Y ); + add( U2, OneAndHalf, Z ); + sub( Y, X, X ); + div( YY1, OneAndHalf, T ); + div( YY1, Y, Y ); + add( Z, U2, t ); + sub( t, T, T ); + sub( Z, Y, Y ); + div( Y2, Z, Z ); + add( Y2, U2, YY1 ); + div( Y2, YY1, YY1 ); + sub( OneAndHalf, Z, Z ); + sub( Y2, YY1, Y2 ); + sub( U1, F9, YY1 ); + div( F9, YY1, YY1 ); + k = cmp( X, Zero ); + k |= cmp( Y, Zero ); + k |= cmp( Z, Zero ); + k |= cmp( T, Zero ); + k |= cmp( Y2, Zero ); + sub( Half, YY1, t ); + sub( Half, F9, t2 ); + k |= cmp( t, t2 ); + if( k == 0 ) + { + RDiv = Rounded; + printf("Division appears to round correctly.\n"); + if(GDiv == No) + printf("Division test inconsistent\n"); + } + else + { + k = 0; + if( cmp( X, Zero ) >= 0 ) + k = 1; + if( cmp( Y, Zero ) >= 0 ) + k = 1; + if( cmp( Z, Zero ) >= 0 ) + k = 1; + if( cmp( T, Zero ) >= 0 ) + k = 1; + if( cmp( Y, Zero ) >= 0 ) + k = 1; + sub( Half, YY1, t ); + sub( Half, F9, t2 ); + if( cmp( t, t2 ) >= 0 ) + k = 1; + if( k == 0 ) + { + RDiv = Chopped; + printf("Division appears to chop.\n"); + } + } + } +if(RDiv == Other) + printf("/ is neither chopped nor correctly rounded.\n"); +div( Radix, One, BInvrse ); +mul( BInvrse, Radix, t ); +sub( Half, t, t ); +if( cmp( t, Half ) != 0 ) + { + ErrCnt[Failure] += 1; + printf( "Radix * ( 1 / Radix ) differs from 1\n" ); + } + +Milestone = 50; +/*=============================================*/ +add( F9, U1, t ); +sub( Half, t, t ); +k = cmp( t, Half ); +add( BMinusU2, U2, t ); +sub( One, t, t ); +sub( One, Radix, t2 ); +k |= cmp( t, t2 ); +if( k != 0 ) + { + ErrCnt[Failure] += 1; + printf( "Incomplete carry-propagation in Addition\n" ); + } +mul( U1, U1, X ); +sub( X, One, X ); +sub( U2, One, Y ); +mul( U2, Y, Y ); +add( One, Y, Y ); +sub( Half, F9, Z ); +sub( Half, X, X ); +sub( Z, X, X ); +sub( One, Y, Y ); +if( (cmp(X,Zero) == 0) && (cmp(Y,Zero) == 0) ) + { + RAddSub = Chopped; + printf("Add/Subtract appears to be chopped.\n"); + } +if(GAddSub == Yes) + { + add( Half, U2, X ); + mul( X, U2, X ); + sub( U2, Half, Y ); + mul( Y, U2, Y ); + add( One, X, X ); + add( One, Y, Y ); + add( One, U2, t ); + sub( X, t, X ); + sub( Y, One, Y ); + k = cmp(X,Zero); + if( k ) + printf( "1+U2-[u2(1/2+U2)+1] != 0\n" ); + k2 = cmp(Y,Zero); + if( k2 ) + printf( "1-[U2(1/2-U2)+1] != 0\n" ); + k |= k2; + if( k == 0 ) + { + add( Half, U2, X ); + mul( X, U1, X ); + sub( U2, Half, Y ); + mul( Y, U1, Y ); + sub( X, One, X ); + sub( Y, One, Y ); + sub( X, F9, X ); + sub( Y, One, Y ); + k = cmp(X,Zero); + if( k ) + printf( "F9-[1-U1(1/2+U2)] != 0\n" ); + k2 = cmp(Y,Zero); + if( k2 ) + printf( "1-[1-U1(1/2-U2)] != 0\n" ); + k |= k2; + if( k == 0 ) + { + RAddSub = Rounded; + printf("Addition/Subtraction appears to round correctly.\n"); + if(GAddSub == No) + printf( "Add/Subtract test inconsistent\n"); + } + else + { + printf("Addition/Subtraction neither rounds nor chops.\n"); + } + } + else + printf("Addition/Subtraction neither rounds nor chops.\n"); + } +else + printf("Addition/Subtraction neither rounds nor chops.\n"); + +mov( One, S ); +add( One, Half, X ); +mul( Half, X, X ); +add( One, X, X ); +add( One, U2, Y ); +mul( Y, Half, Y ); +sub( Y, X, Z ); +sub( X, Y, T ); +add( Z, T, StickyBit ); +if( cmp(StickyBit, Zero) != 0 ) + { + mov( Zero, S ); + ErrCnt[Flaw] += 1; + printf( "(X - Y) + (Y - X) is non zero!\n" ); + } +mov( Zero, StickyBit ); +FLOOR( RadixD2, t ); +k2 = cmp( t, RadixD2 ); +k = 1; +if( (GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes) + && (RMult == Rounded) && (RDiv == Rounded) + && (RAddSub == Rounded) && (k2 == 0) ) + { + printf("Checking for sticky bit.\n"); + k = 0; + add( Half, U1, X ); + mul( X, U2, X ); + mul( Half, U2, Y ); + add( One, Y, Z ); + add( One, X, T ); + sub( One, Z, t ); + sub( One, T, t2 ); + if( cmp(t,Zero) > 0 ) + { + k = 1; + printf( "[1+(1/2)U2]-1 > 0\n" ); + } + if( cmp(t2,U2) < 0 ) + { + k = 1; + printf( "[1+U2(1/2+U1)]-1 < U2\n" ); + } + add( T, Y, Z ); + sub( X, Z, Y ); + sub( T, Z, t ); + sub( T, Y, t2 ); + if( cmp(t,U2) < 0 ) + { + k = 1; + printf( "[[1+U2(1/2+U1)]+(1/2)U2]-[1+U2(1/2+U1)] < U2\n" ); + } + if( cmp(t2,Zero) != 0 ) + { + k = 1; + printf( "(1/2)U2-[1+U2(1/2+U1)] != 0\n" ); + } + add( Half, U1, X ); + mul( X, U1, X ); + mul( Half, U1, Y ); + sub( Y, One, Z ); + sub( X, One, T ); + sub( One, Z, t ); + sub( F9, T, t2 ); + if( cmp(t,Zero) != 0 ) + { + k = 1; + printf( "(1-(1/2)U1)-1 != 0\n" ); + } + if( cmp(t2,Zero) != 0 ) + { + k = 1; + printf( "[1-U1(1/2+U1)]-F9 != 0\n" ); + } + sub( U1, Half, Z ); + mul( Z, U1, Z ); + sub( Z, F9, T ); + sub( Y, F9, Q ); + sub( F9, T, t ); + if( cmp( t, Zero ) != 0 ) + { + k = 1; + printf( "[F9-U1(1/2-U1)]-F9 != 0\n" ); + } + sub( U1, F9, t ); + sub( Q, t, t ); + if( cmp( t, Zero ) != 0 ) + { + k = 1; + printf( "(F9-U1)-(F9-(1/2)U1) != 0\n" ); + } + add( One, U2, Z ); + mul( Z, OneAndHalf, Z ); + add( OneAndHalf, U2, T ); + sub( Z, T, T ); + add( U2, T, T ); + div( Radix, Half, X ); + add( One, X, X ); + mul( Radix, U2, Y ); + add( One, Y, Y ); + mul( X, Y, Z ); + if( cmp( T, Zero ) != 0 ) + { + k = 1; + printf( "(3/2+U2)-3/2(1+U2)+U2 != 0\n" ); + } + mul( Radix, U2, t ); + add( X, t, t ); + sub( Z, t, t ); + if( cmp( t, Zero ) != 0 ) + { + k = 1; + printf( "(1+1/2Radix)+Radix*U2-[1+1/(2Radix)][1+Radix*U2] != 0\n" ); + } + if( cmp(Radix, Two) != 0 ) + { + add( Two, U2, X ); + div( Two, X, Y ); + sub( One, Y, t ); + if( cmp( t, Zero) != 0 ) + k = 1; + } + } +if( k == 0 ) + { + printf("Sticky bit apparently used correctly.\n"); + mov( One, StickyBit ); + } +else + { + printf("Sticky bit used incorrectly or not at all.\n"); + } + +if( GMult == No || GDiv == No || GAddSub == No || + RMult == Other || RDiv == Other || RAddSub == Other) + { + ErrCnt[Flaw] += 1; + printf("lack(s) of guard digits or failure(s) to correctly round or chop\n"); +printf( "(noted above) count as one flaw in the final tally below\n" ); + } +/*=============================================*/ +Milestone = 60; +/*=============================================*/ +printf("\n"); +printf("Does Multiplication commute? "); +printf("Testing on %d random pairs.\n", NoTrials); +SQRT( Three, Random9 ); +mov( Third, Random1 ); +I = 1; +do + { + Random(); + mov( Random1, X ); + Random(); + mov( Random1, Y ); + mul( Y, X, Z9 ); + mul( X, Y, Z ); + sub( Z9, Z, Z9 ); + I = I + 1; + } +while ( ! ((I > NoTrials) || (cmp(Z9,Zero) != 0))); +if(I == NoTrials) + { + div( Three, Half, t ); + add( One, t, Random1 ); + add( U2, U1, t ); + add( t, One, Random2 ); + mul( Random1, Random2, Z ); + mul( Random2, Random1, Y ); +/* Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half / + * Three) * ((U2 + U1) + One); + */ + div( Three, Half, t2 ); + add( One, t2, t2 ); + add( U2, U1, t ); + add( t, One, t ); + mul( t2, t, Z9 ); + mul( t2, t, t ); + sub( t, Z9, Z9 ); + } +if(! ((I == NoTrials) || (cmp(Z9,Zero) == 0))) + { + ErrCnt[Defect] += 1; + printf( "X * Y == Y * X trial fails.\n"); + } +else + { + printf(" No failures found in %d integer pairs.\n", NoTrials); + } +/*=============================================*/ +Milestone = 70; +/*=============================================*/ +sqtest(); +Milestone = 90; +pow1test(); + +Milestone = 110; + +printf("Seeking Underflow thresholds UfThold and E0.\n"); +mov( U1, D ); +FLOOR( Precision, t ); +if( cmp(Precision, t) != 0 ) + { + mov( BInvrse, D ); + mov( Precision, X ); + do + { + mul( D, BInvrse, D ); + sub( One, X, X ); + } + while( cmp(X, Zero) > 0 ); + } +mov( One, Y ); +mov( D, Z ); +/* ... D is power of 1/Radix < 1. */ +sigsave = sigfpe; +if( setjmp(ovfl_buf) ) + goto under0; +do + { + mov( Y, C ); + mov( Z, Y ); + mul( Y, Y, Z ); + add( Z, Z, t ); + } +while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) ); + +under0: +sigsave = 0; + +mov( C, Y ); +mul( Y, D, Z ); +sigsave = sigfpe; +if( setjmp(ovfl_buf) ) + goto under1; +do + { + mov( Y, C ); + mov( Z, Y ); + mul( Y, D, Z ); + add( Z, Z, t ); + } +while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) ); + +under1: +sigsave = 0; + +if( cmp(Radix,Two) < 0 ) + mov( Two, HInvrse ); +else + mov( Radix, HInvrse ); +div( HInvrse, One, H ); +/* ... 1/HInvrse == H == Min(1/Radix, 1/2) */ +div( C, One, CInvrse ); +mov( C, E0 ); +mul( E0, H, Z ); +/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */ +sigsave = sigfpe; +if( setjmp(ovfl_buf) ) + goto under2; +do + { + mov( E0, Y ); + mov( Z, E0 ); + mul( E0, H, Z ); + add( Z, Z, t ); + } +while( (cmp(E0,Z) > 0) && (cmp(t,Z) > 0) ); + +under2: +sigsave = 0; + +mov( E0, UfThold ); +mov( Zero, E1 ); +mov( Zero, Q ); +mov( U2, E9 ); +add( One, E9, S ); +mul( C, S, D ); +if( cmp(D,C) <= 0 ) + { + mul( Radix, U2, E9 ); + add( One, E9, S ); + mul( C, S, D ); + if( cmp(D, C) <= 0 ) + { + ErrCnt[Failure] += 1; + printf( "multiplication gets too many last digits wrong.\n" ); + mov( E0, Underflow ); + mov( Zero, YY1 ); + mov( Z, PseudoZero ); + } + } +else + { + mov( D, Underflow ); + mul( Underflow, H, PseudoZero ); + mov( Zero, UfThold ); + do + { + mov( Underflow, YY1 ); + mov( PseudoZero, Underflow ); + add( E1, E1, t ); + if( cmp(t, E1) <= 0) + { + mul( Underflow, HInvrse, Y2 ); + sub( Y2, YY1, E1 ); + FABS( E1 ); + mov( YY1, Q ); + if( (cmp( UfThold, Zero ) == 0) + && (cmp(YY1, Y2) != 0) ) + mov( YY1, UfThold ); + } + mul( PseudoZero, H, PseudoZero ); + add( PseudoZero, PseudoZero, t ); + } + while( (cmp(Underflow, PseudoZero) > 0) + && (cmp(t, PseudoZero) > 0) ); + } +/* Comment line 4530 .. 4560 */ +if( cmp(PseudoZero, Zero) != 0 ) + { + printf("\n"); + mov(PseudoZero, Z ); +/* ... Test PseudoZero for "phoney- zero" violates */ +/* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero + ... */ + if( cmp(PseudoZero, Zero) <= 0 ) + { + ErrCnt[Failure] += 1; + printf("Positive expressions can underflow to an\n"); + printf("allegedly negative value\n"); + printf("PseudoZero that prints out as: " ); + show( PseudoZero ); + mov( PseudoZero, X ); + neg( X ); + if( cmp(X, Zero) <= 0 ) + { + printf("But -PseudoZero, which should be\n"); + printf("positive, isn't; it prints out as " ); + show( X ); + } + } + else + { + ErrCnt[Flaw] += 1; + printf( "Underflow can stick at an allegedly positive\n"); + printf("value PseudoZero that prints out as " ); + show( PseudoZero ); + } +/* TstPtUf();*/ + } + +/*=============================================*/ +Milestone = 120; +/*=============================================*/ +mul( CInvrse, Y, t ); +mul( CInvrse, YY1, t2 ); +if( cmp(t,t2) > 0 ) + { + mul( H, S, S ); + mov( Underflow, E0 ); + } +if(! ((cmp(E1,Zero) == 0) || (cmp(E1,E0) == 0)) ) + { + ErrCnt[Defect] += 1; + if( cmp(E1,E0) < 0 ) + { + printf("Products underflow at a higher"); + printf(" threshold than differences.\n"); + if( cmp(PseudoZero,Zero) == 0 ) + mov( E1, E0 ); + } + else + { + printf("Difference underflows at a higher"); + printf(" threshold than products.\n"); + } + } +printf("Smallest strictly positive number found is E0 = " ); +show( E0 ); +mov( E0, Z ); +TstPtUf(); +mov( E0, Underflow ); +if(N == 1) + mov( Y, Underflow ); +I = 4; +if( cmp(E1,Zero) == 0 ) + I = 3; +if( cmp( UfThold,Zero) == 0 ) + I = I - 2; +UfNGrad = True; +switch(I) + { + case 1: + mov( Underflow, UfThold ); + mul( CInvrse, Q, t ); + mul( CInvrse, Y, t2 ); + mul( t2, S, t2 ); + if( cmp( t, t2 ) != 0 ) + { + mov( Y, UfThold ); + ErrCnt[Failure] += 1; + printf( "Either accuracy deteriorates as numbers\n"); + printf("approach a threshold = " ); + show( UfThold ); + printf(" coming down from " ); + show( C ); + printf(" or else multiplication gets too many last digits wrong.\n"); + } + break; + + case 2: + ErrCnt[Failure] += 1; + printf( "Underflow confuses Comparison which alleges that\n"); + printf("Q == Y while denying that |Q - Y| == 0; these values\n"); + printf("print out as Q = " ); + show( Q ); + printf( ", Y = " ); + show( Y ); + sub( Y2, Q, t ); + FABS(t); + printf ("|Q - Y| = " ); + show( t ); + mov( Q, UfThold ); + break; + + case 3: + mov( X, X ); + break; + + case 4: + div( E9, E1, t ); + sub( t, UfThold, t ); + FABS(t); + if( (cmp(Q,UfThold) == 0) && (cmp(E1,E0) == 0) + && (cmp(t,E1) <= 0) ) + { + UfNGrad = False; + printf("Underflow is gradual; it incurs Absolute Error =\n"); + printf("(roundoff in UfThold) < E0.\n"); + mul( E0, CInvrse, Y ); + add( OneAndHalf, U2, t ); + mul( Y, t, Y ); + add( One, U2, X ); + mul( CInvrse, X, X ); + div( X, Y, t ); + IEEE = (cmp(t,E0) == 0); + if( IEEE == 0 ) + { + printf( "((CInvrse E0) (1.5+U2)) / (CInvrse (1+U2)) != E0\n" ); + printf( "CInvrse = " ); + show( CInvrse ); + printf( "E0 = " ); + show( E0 ); + printf( "U2 = " ); + show( U2 ); + printf( "X = " ); + show(X); + printf( "Y = " ); + show(Y); + printf( "Y/X = " ); + show(t); + } + } + } +if(UfNGrad) + { + printf("\n"); + div( UfThold, Underflow, R ); + SQRT( R, R ); + if( cmp(R,H) <= 0) + { + mul( R, UfThold, Z ); +/* X = Z * (One + R * H * (One + H));*/ + add( One, H, X ); + mul( H, X, X ); + mul( R, X, X ); + add( One, X, X ); + mul( Z, X, X ); + } + else + { + mov( UfThold, Z ); +/*X = Z * (One + H * H * (One + H));*/ + add( One, H, X ); + mul( H, X, X ); + mul( H, X, X ); + add( One, X, X ); + mul( Z, X, X ); + } + sub( Z, X, t ); +/* if(! ((cmp(X,Z) == 0) || (cmp(t,Zero) != 0)) )*/ + if( (cmp(X,Z) != 0) && (cmp(t,Zero) == 0) ) + { +/* ErrCnt[Flaw] += 1;*/ + ErrCnt[Serious] += 1; + printf("X = " ); + show( X ); + printf( "\tis not equal to Z = " ); + show( Z ); +/* sub( Z, X, Z9 );*/ + printf("yet X - Z yields " ); + show( t ); + printf("which compares equal to " ); + show( Zero ); + printf(" Should this NOT signal Underflow, "); + printf("this is a SERIOUS DEFECT\nthat causes "); + printf("confusion when innocent statements like\n");; + printf(" if (X == Z) ... else"); + printf(" ... (f(X) - f(Z)) / (X - Z) ...\n"); + printf("encounter Division by Zero although actually\n"); + printf("X / Z = 1 + " ); + div( Z, X, t ); + sub( Half, t, t ); + sub( Half, t, t ); + show(t); + } + } +printf("The Underflow threshold is " ); +show( UfThold ); +printf( "below which calculation may suffer larger Relative error than" ); +printf( " merely roundoff.\n"); +mul( U1, U1, Y2 ); +mul( Y2, Y2, Y ); +mul( Y, U1, Y2 ); +if( cmp( Y2,UfThold) <= 0 ) + { + if( cmp(Y,E0) > 0 ) + { + ErrCnt[Defect] += 1; + I = 5; + } + else + { + ErrCnt[Serious] += 1; + I = 4; + } + printf("Range is too narrow; U1^%d Underflows.\n", I); + } +Milestone = 130; + +/*Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;*/ +LOG( UfThold, Y ); +LOG( HInvrse, t ); +div( t, Y, Y ); +mul( TwoForty, Y, Y ); +sub( Y, Half, Y ); +FLOOR( Y, Y ); +div( TwoForty, Y, Y ); +neg(Y); +sub( One, Y, Y2 ); /* ***** changed from Y2 = Y + Y */ +printf("Since underflow occurs below the threshold\n"); +printf("UfThold = " ); +show( HInvrse ); +printf( "\tto the power " ); +show( Y ); +printf( "only underflow should afflict the expression " ); +show( HInvrse ); +printf( "\tto the power " ); +show( Y2 ); +POW( HInvrse, Y2, V9 ); +printf("Actually calculating yields: " ); +show( V9 ); +add( Radix, Radix, t ); +add( t, E9, t ); +mul( t, UfThold, t ); +if( (cmp(V9,Zero) < 0) || (cmp(V9,t) > 0) ) + { + ErrCnt[Serious] += 1; + printf( "this is not between 0 and underflow\n"); + printf(" threshold = " ); + show( UfThold ); + } +else + { + add( One, E9, t ); + mul( UfThold, t, t ); + if( cmp(V9,t) <= 0 ) + printf("This computed value is O.K.\n"); + else + { + ErrCnt[Defect] += 1; + printf( "this is not between 0 and underflow\n"); + printf(" threshold = " ); + show( UfThold ); + } + } + +Milestone = 140; + +pow2test(); + +/*=============================================*/ +Milestone = 160; +/*=============================================*/ +Pause(); +printf("Searching for Overflow threshold:\n"); +printf("This may generate an error.\n"); +sigsave = sigfpe; +I = 0; +mov( CInvrse, Y ); /* a large power of 2 */ +neg(Y); +mul( HInvrse, Y, V9 ); /* HInvrse = 2 */ +if (setjmp(ovfl_buf)) + goto overflow; +do + { + mov( Y, V ); + mov( V9, Y ); + mul( HInvrse, Y, V9 ); + } +while( cmp(V9,Y) < 0 ); /* V9 = 2 * Y */ +I = 1; + +overflow: + +show( HInvrse ); +printf( "\ttimes " ); +show( Y ); +printf( "\tequals " ); +show( V9 ); + +mov( V9, Z ); +printf("Can `Z = -Y' overflow?\n"); +printf("Trying it on Y = " ); +show(Y); +mov( Y, V9 ); +neg( V9 ); +mov( V9, V0 ); +sub( Y, V, t ); +add( V, V0, t2 ); +if( cmp(t,t2) == 0 ) + printf("Seems O.K.\n"); +else + { + printf("finds a Flaw, -(-Y) differs from Y.\n"); + printf( "V-Y=t:" ); + show(V); + show(Y); + show(t); + printf( "V+V0=t2:" ); + show(V); + show(V0); + show(t2); + ErrCnt[Flaw] += 1; + } +if( (cmp(Z, Y) != 0) && (I != 0) ) + { + ErrCnt[Serious] += 1; + printf("overflow past " ); + show( Y ); + printf( "\tshrinks to " ); + show( Z ); + printf( "= Y * " ); + show( HInvrse ); + } +/*Y = V * (HInvrse * U2 - HInvrse);*/ +mul( HInvrse, U2, Y ); +sub( HInvrse, Y, Y ); +mul( V, Y, Y ); +/*Z = Y + ((One - HInvrse) * U2) * V;*/ +sub( HInvrse, One, Z ); +mul( Z, U2, Z ); +mul( Z, V, Z ); +add( Y, Z, Z ); +if( cmp(Z,V0) < 0 ) + mov( Z, Y ); +if( cmp(Y,V0) < 0) + mov( Y, V ); +sub( V, V0, t ); +if( cmp(t,V0) < 0 ) + mov( V0, V ); +printf("Overflow threshold is V = " ); +show( V ); +if(I) + { + printf("Overflow saturates at V0 = " ); + show( V0 ); + } +else +printf("There is no saturation value because the system traps on overflow.\n"); + +mul( V, One, V9 ); +printf("No Overflow should be signaled for V * 1 = " ); +show( V9 ); +div( One, V, V9 ); + printf(" nor for V / 1 = " ); + show( V9 ); + printf("Any overflow signal separating this * from the one\n"); + printf("above is a DEFECT.\n"); +/*=============================================*/ +Milestone = 170; +/*=============================================*/ +mov( V, t ); +neg( t ); +k = 0; +if( cmp(t,V) >= 0 ) + k = 1; +mov( V0, t ); +neg( t ); +if( cmp(t,V0) >= 0 ) + k = 1; +mov( UfThold, t ); +neg(t); +if( cmp(t,V) >= 0 ) + k = 1; +if( cmp(UfThold,V) >= 0 ) + k = 1; +if( k != 0 ) + { + ErrCnt[Failure] += 1; + printf( "Comparisons involving +-"); + show( V ); + show( V0 ); + show( UfThold ); + printf("are confused by Overflow." ); + } +/*=============================================*/ +Milestone = 175; +/*=============================================*/ +printf("\n"); +for(Indx = 1; Indx <= 3; ++Indx) { + switch(Indx) + { + case 1: mov(UfThold, Z); break; + case 2: mov( E0, Z); break; + case 3: mov(PseudoZero, Z); break; + } +if( cmp(Z, Zero) != 0 ) + { + SQRT( Z, V9 ); + mul( V9, V9, Y ); + mul( Radix, E9, t ); + sub( t, One, t ); + div( t, Y, t ); + add( One, Radix, t2 ); + add( t2, E9, t2 ); + mul( t2, Z, t2 ); + if( (cmp(t,Z) < 0) || (cmp(Y,t2) > 0) ) + { + if( cmp(V9,U1) > 0 ) + ErrCnt[Serious] += 1; + else + ErrCnt[Defect] += 1; + printf("Comparison alleges that what prints as Z = " ); + show( Z ); + printf(" is too far from sqrt(Z) ^ 2 = " ); + show( Y ); + } + } +} + +Milestone = 180; + +for(Indx = 1; Indx <= 2; ++Indx) + { + if(Indx == 1) + mov( V, Z ); + else + mov( V0, Z ); + SQRT( Z, V9 ); + mul( Radix, E9, X ); + sub( X, One, X ); + mul( X, V9, X ); + mul( V9, X, V9 ); + mul( Two, Radix, t ); + mul( t, E9, t ); + sub( t, One, t ); + mul( t, Z, t ); + if( (cmp(V9,t) < 0) || (cmp(V9,Z) > 0) ) + { + mov( V9, Y ); + if( cmp(X,W) < 0 ) + ErrCnt[Serious] += 1; + else + ErrCnt[Defect] += 1; + printf("Comparison alleges that Z = " ); + show( Z ); + printf(" is too far from sqrt(Z) ^ 2 :" ); + show( Y ); + } + } + +Milestone = 190; + +Pause(); +mul( UfThold, V, X ); +mul( Radix, Radix, Y ); +mul( X, Y, t ); +if( (cmp(t,One) < 0) || (cmp(X,Y) > 0) ) + { + mul( X, Y, t ); + div( U1, Y, t2 ); + if( (cmp(t,U1) < 0) || (cmp(X,t2) > 0) ) + { + ErrCnt[Defect] += 1; + printf( "Badly " ); + } + else + { + ErrCnt[Flaw] += 1; + } + printf(" unbalanced range; UfThold * V = " ); + show( X ); + printf( "\tis too far from 1.\n"); + } +Milestone = 200; + +for(Indx = 1; Ind