Add in some math lib tests

1 files changed, 3550 insertions, 0 deletions

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;

+

+/* 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; Indx <= 5; ++Indx)

+	{

+	mov( F9, X );

+	switch(Indx)

+		{

+		case 2: add( One, U2, X ); break;

+		case 3: mov( V, X ); break;

+		case 4: mov(UfThold,X); break;

+		case 5: mov(Radix,X);

+		}

+	mov( X, Y );

+

+	sigsave = sigfpe;

+	if (setjmp(ovfl_buf))

+		{

+		printf("  X / X  traps when X = " );

+		show( X );

+		}

+	else

+		{

+/*V9 = (Y / X - Half) - Half;*/

+		div( X, Y, t );

+		sub( Half, t, t );

+		sub( Half, t, V9 );

+		if( cmp(V9,Zero) == 0 )

+			continue;

+		mov( U1, t );

+		neg(t);

+		if( (cmp(V9,t) == 0) && (Indx < 5) )

+			{

+			ErrCnt[Flaw] += 1;

+			}

+		else

+			{

+			ErrCnt[Serious] += 1;

+			}

+		printf("  X / X differs from 1 when X = " );

+		show( X );

+		printf("  instead, X / X - 1/2 - 1/2 = " );

+		show( V9 );

+		}

+	}

+

+	Pause();

+	printf("\n");

+	{

+		static char *msg[] = {

+			"FAILUREs  encountered =",

+			"SERIOUS DEFECTs  discovered =",

+			"DEFECTs  discovered =",

+			"FLAWs  discovered =" };

+		int i;

+		for(i = 0; i < 4; i++) if (ErrCnt[i])

+			printf("The number of  %-29s %d.\n",

+				msg[i], ErrCnt[i]);

+		}

+	printf("\n");

+	if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]

+			+ ErrCnt[Flaw]) > 0) {

+		if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[

+			Defect] == 0) && (ErrCnt[Flaw] > 0)) {

+			printf("The arithmetic diagnosed seems ");

+			printf("satisfactory though flawed.\n");

+			}

+		if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)

+			&& ( ErrCnt[Defect] > 0)) {

+			printf("The arithmetic diagnosed may be acceptable\n");

+			printf("despite inconvenient Defects.\n");

+			}

+		if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {

+			printf("The arithmetic diagnosed has ");

+			printf("unacceptable serious defects.\n");

+			}

+		if (ErrCnt[Failure] > 0) {

+			printf("Fatal FAILURE may have spoiled this");

+			printf(" program's subsequent diagnoses.\n");

+			}

+		}

+	else {

+		printf("No failures, defects nor flaws have been discovered.\n");

+		if (! ((RMult == Rounded) && (RDiv == Rounded)

+			&& (RAddSub == Rounded) && (RSqrt == Rounded))) 

+			printf("The arithmetic diagnosed seems satisfactory.\n");

+		else {

+			k = 0;

+			if( cmp( Radix, Two ) == 0 )

+				k = 1;

+			if( cmp( Radix, Ten ) == 0 )

+				k = 1;

+			if( (cmp(StickyBit,One) >= 0) && (k == 1) )

+				{

+				printf("Rounding appears to conform to ");

+				printf("the proposed IEEE standard P");

+				k = 0;

+				k |= cmp( Radix, Two );

+				mul( Four, Three, t );

+				mul( t, Two, t );

+				sub( t, Precision, t );

+				sub( TwentySeven, Precision, t2 );

+				sub( TwentySeven, t2, t2 );

+				add( t2, One, t2 );

+				mul( t2, t, t );

+				if( (cmp(Radix,Two) == 0)

+					&& (cmp(t,Zero) == 0) )

+					printf("754");

+				else

+					printf("854");

+				if(IEEE)

+					printf(".\n");

+				else

+					{

+			printf(",\nexcept for possibly Double Rounding");

+			printf(" during Gradual Underflow.\n");

+					}

+				}

+		printf("The arithmetic diagnosed appears to be excellent!\n");

+			}

+		}

+	if (fpecount)

+		printf("\nA total of %d floating point exceptions were registered.\n",

+			fpecount);

+	printf("END OF TEST.\n");

+	}

+

+

+/* Random */

+/*  Random computes

+     X = (Random1 + Random9)^5

+     Random1 = X - FLOOR(X) + 0.000005 * X;

+   and returns the new value of Random1

+*/

+

+

+static int randflg = 0;

+FLOAT(C5em6);

+

+Random()

+{

+

+if( randflg == 0 )

+	{

+	mov( Six, t );

+	neg(t);

+	POW( Ten, t, t );

+	mul( Five, t, C5em6 );

+	randflg = 1;

+	}

+add( Random1, Random9, t );

+mul( t, t, t2 );

+mul( t2, t2, t2 );

+mul( t, t2, t );

+FLOOR(t, t2 );

+sub( t2, t, t2 );

+mul( t, C5em6, t );

+add( t, t2, Random1 );

+/*return(Random1);*/

+}

+

+/* SqXMinX */

+

+SqXMinX( ErrKind )

+int ErrKind;

+{

+mul( X, BInvrse, t2 );

+sub( t2, X, t );

+/*SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;*/

+mul( X, X, Sqarg );

+SQRT( Sqarg, SqEr );

+sub( t2, SqEr, SqEr );

+sub( t, SqEr, SqEr );

+div( OneUlp, SqEr, SqEr );

+if( cmp(SqEr,Zero) != 0)

+	{

+	Showsq( 0 );

+	add( J, One, J );

+	ErrCnt[ErrKind] += 1;

+	printf("sqrt of " );

+	mul( X, X, t );

+	show( t );

+	printf( "minus " );

+	show( X );

+	printf( "equals " );

+	mul( OneUlp, SqEr, t );

+	show( t );

+	printf("\tinstead of correct value 0 .\n");

+	}

+}

+

+/* NewD */

+

+NewD()

+{

+mul( Z1, Q, X );

+/*X = FLOOR(Half - X / Radix) * Radix + X;*/

+div( Radix, X, t );

+sub( t, Half, t );

+FLOOR( t, t );

+mul( t, Radix, t );

+add( t, X, X );

+/*Q = (Q - X * Z) / Radix + X * X * (D / Radix);*/

+mul( X, Z, t );

+sub( t, Q, t );

+div( Radix, t, t );

+div( Radix, D, t2 );

+mul( X, t2, t2 );

+mul( X, t2, t2 );

+add( t, t2, Q );

+/*Z = Z - Two * X * D;*/

+mul( Two, X, t );

+mul( t, D, t );

+sub( t, Z, Z );

+

+if( cmp(Z,Zero) <= 0)

+	{

+	neg(Z);

+	neg(Z1);

+	}

+mul( Radix, D, D );

+}

+

+/* SR3750 */

+

+SR3750()

+{

+sub( Radix, X, t );

+sub( Radix, Z2, t2 );

+k = 0;

+if( cmp(t,t2) < 0 )

+	k = 1;

+sub( Z2, X, t );

+sub( Z2, W, t2 );

+if( cmp(t,t2) > 0 )

+	k = 1;

+/*if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {*/

+if( k == 0 )

+	{

+	I = I + 1;

+	mul( X, D, X2 );

+	mov( X2, Sqarg );

+	SQRT( X2, X2 );

+/*Y2 = (X2 - Z2) - (Y - Z2);*/

+	sub( Z2, X2, Y2 );

+	sub( Z2, Y, t );

+	sub( t, Y2, Y2 );

+	sub( Half, Y, X2 );

+	div( X2, X8, X2 );

+	mul( Half, X2, t );

+	mul( t, X2, t );

+	sub( t, X2, X2 );

+/*SqEr = (Y2 + Half) + (Half - X2);*/

+	add( Y2, Half, SqEr );

+	sub( X2, Half, t );

+	add( t, SqEr, SqEr );

+	Showsq( -1 );

+	sub( X2, Y2, SqEr );

+	Showsq( 1 );

+	}

+}

+

+/* IsYeqX */

+

+IsYeqX()

+{

+if( cmp(Y,X) != 0 )

+	{

+	if (N <= 0)

+		{

+		if( (cmp(Z,Zero) == 0) && (cmp(Q,Zero) <= 0) )

+			printf("WARNING:  computing\n");

+		else

+			{

+			ErrCnt[Defect] += 1;

+			printf( "computing\n");

+			}

+		show( Z );

+		printf( "\tto the power " );

+		show( Q );

+		printf("\tyielded " );

+		show( Y );

+		printf("\twhich compared unequal to correct " );

+		show( X );

+		sub( X, Y, t );

+		printf("\t\tthey differ by " );

+		show( t );

+		}

+	N = N + 1; /* ... count discrepancies. */

+	}

+}

+

+/* SR3980 */

+

+SR3980()

+{

+long li;

+

+do

+	{

+/*Q = (FLOAT) I;*/

+	li = I;

+	LTOF( &li, Q );

+	POW( Z, Q, Y );

+	IsYeqX();

+	if(++I > M)

+		break;

+	mul( Z, X, X );

+	}

+while( cmp(X,W) < 0 );

+}

+

+/* PrintIfNPositive */

+

+PrintIfNPositive()

+{

+if(N > 0)

+	printf("Similar discrepancies have occurred %d times.\n", N);

+}

+

+

+/* TstPtUf */

+

+TstPtUf()

+{

+N = 0;

+if( cmp(Z,Zero) != 0)

+	{

+	printf( "Z = " );

+	show(Z);

+	printf("Since comparison denies Z = 0, evaluating ");

+	printf("(Z + Z) / Z should be safe.\n");

+	sigsave = sigfpe;

+	if (setjmp(ovfl_buf))

+		goto very_serious;

+	add( Z, Z, Q9 );

+	div( Z, Q9, Q9 );

+	printf("What the machine gets for (Z + Z) / Z is " );

+	show( Q9 );

+	sub( Two, Q9, t );

+	FABS(t);

+	mul( Radix, U2, t2 );

+	if( cmp(t,t2) < 0 )

+		{

+		printf("This is O.K., provided Over/Underflow");

+		printf(" has NOT just been signaled.\n");

+		}

+	else

+		{

+		if( (cmp(Q9,One) < 0) || (cmp(Q9,Two) > 0) )

+			{

+very_serious:

+			N = 1;

+			ErrCnt [Serious] = ErrCnt [Serious] + 1;

+			printf("This is a VERY SERIOUS DEFECT!\n");

+			}

+		else

+			{

+			N = 1;

+			ErrCnt[Defect] += 1;

+			printf("This is a DEFECT!\n");

+			}

+		}

+	mul( Z, One, V9 );

+	mov( V9, Random1 );

+	mul( One, Z, V9 );

+	mov( V9, Random2 );

+	div( One, Z, V9 );

+	if( (cmp(Z,Random1) == 0) && (cmp(Z,Random2) == 0)

+		&& (cmp(Z,V9) == 0) )

+		{

+		if (N > 0)

+			Pause();

+		}

+	else

+		{

+		N = 1;

+		ErrCnt[Defect] += 1;

+		printf( "What prints as Z = ");

+		show( Z );

+		printf( "\tcompares different from " );

+		if( cmp(Z,Random1) != 0)

+			{

+			printf("Z * 1 = " );

+			show( Random1 );

+			}

+		if( (cmp(Z,Random2) != 0)

+			|| (cmp(Random2,Random1) != 0) )

+			{

+			printf("1 * Z == " );

+			show( Random2 );

+			}

+		if( cmp(Z,V9) != 0 )

+			{

+			printf("Z / 1 = " );

+			show( V9 );

+			}

+		if( cmp(Random2,Random1) != 0 )

+			{

+			ErrCnt[Defect] += 1;

+			printf( "Multiplication does not commute!\n");

+			printf("\tComparison alleges that 1 * Z = " );

+			show(Random2);

+			printf("\tdiffers from Z * 1 = " );

+			show(Random1);

+			}

+		Pause();

+		}

+	}

+}

+

+Pause()

+{

+}

+

+Sign( x, y )

+FSIZE *x, *y;

+{

+

+if( cmp( x, Zero ) < 0 )

+	{

+	mov( One, y );

+	neg( y );

+	}

+else

+	{

+	mov( One, y );

+	}

+}

+

+sqtest()

+{

+printf("\nRunning test of square root(x).\n");

+

+RSqrt = Other;

+k = 0;

+SQRT( Zero, t );

+k |= cmp( Zero, t );

+mov( Zero, t );

+neg(t);

+SQRT( t, t2 );

+k |= cmp( t, t2 );

+SQRT( One, t );

+k |= cmp( One, t );

+if( k != 0 )

+ 	{

+	ErrCnt[Failure] += 1;

+	printf( "Square root of 0.0, -0.0 or 1.0 wrong\n");

+	}

+mov( Zero, MinSqEr );

+mov( Zero, MaxSqEr );

+mov( Zero, J );

+mov( Radix, X );

+mov( U2, OneUlp );

+SqXMinX( Serious );

+mov( BInvrse, X );

+mul( BInvrse, U1, OneUlp );

+SqXMinX( Serious );

+mov( U1, X );

+mul( U1, U1, OneUlp );

+SqXMinX( Serious );

+if( cmp(J,Zero) != 0)

+	Pause();

+printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);

+mov( Zero, J );

+mov( Two, X );

+mov( Radix, Y );

+if( cmp(Radix,One) != 0 )

+	{

+	lngint = NoTrials;

+	LTOF( &lngint, t );

+	FTOL( t, &lng2, X );

+	if( lngint != lng2 )

+		{

+		printf( "Integer conversion error\n" );

+		exit(1);

+		}

+	do

+		{

+		mov( Y, X );

+		mul( Radix, Y, Y );

+		sub( X, Y, t2 );

+		}

+	while( ! (cmp(t2,t) >= 0) );

+	}

+mul( X, U2, OneUlp );

+I = 1;

+while(I < 10)

+	{

+	add( X, One, X );

+	SqXMinX( Defect );

+	if( cmp(J,Zero) > 0 )

+		break;

+	I = I + 1;

+	}

+printf("Test for sqrt monotonicity.\n");

+I = - 1;

+mov( BMinusU2, X );

+mov( Radix, Y );

+mul( Radix, U2, Z );

+add( Radix, Z, Z );

+NotMonot = False;

+Monot = False;

+while( ! (NotMonot || Monot))

+	{

+	I = I + 1;

+	SQRT(X, X);

+	SQRT(Y,Q);

+	SQRT(Z,Z);

+	if( (cmp(X,Q) > 0) || (cmp(Q,Z) > 0) )

+		NotMonot = True;

+	else

+		{

+		add( Q, Half, Q );

+		FLOOR( Q, Q );

+		mul( Q, Q, t );

+		if( (I > 0) || (cmp(Radix,t) == 0) )

+			Monot = True;

+		else if (I > 0)

+			{

+			if(I > 1)

+				Monot = True;

+			else

+				{

+				mul( Y, BInvrse, Y );

+				sub( U1, Y, X );

+				add( Y, U1, Z );

+				}

+			}

+		else

+			{

+			mov( Q, Y );

+			sub( U2, Y, X );

+			add( Y, U2, Z );

+			}

+		}

+	}

+if( Monot )

+	printf("sqrt has passed a test for Monotonicity.\n");

+else

+	{

+	ErrCnt[Defect] += 1;

+	printf("sqrt(X) is non-monotonic for X near " );

+	show(Y);

+	}

+/*=============================================*/

+Milestone = 80;

+/*=============================================*/

+add( MinSqEr, Half, MinSqEr );

+sub( Half, MaxSqEr, MaxSqEr);

+/*Y = (SQRT(One + U2) - One) / U2;*/

+add( One, U2, Sqarg );

+SQRT( Sqarg, Y );

+sub( One, Y, Y );

+div( U2, Y, Y );

+/*SqEr = (Y - One) + U2 / Eight;*/

+sub( One, Y, t );

+div( Eight, U2, SqEr );

+add( t, SqEr, SqEr );

+Showsq( 1 );

+div( Eight, U2, SqEr );

+add( Y, SqEr, SqEr );

+Showsq( -1 );

+/*Y = ((SQRT(F9) - U2) - (One - U2)) / U1;*/

+mov( F9, Sqarg );

+SQRT( Sqarg, Y );

+sub( U2, Y, Y );

+sub( U2, One, t );

+sub( t, Y, Y );

+div( U1, Y, Y );

+div( Eight, U1, SqEr );

+add( Y, SqEr, SqEr );

+Showsq( 1 );

+/*SqEr = (Y + One) + U1 / Eight;*/

+div( Eight, U1, t );

+add( Y, One, SqEr );

+add( SqEr, t, SqEr );

+Showsq( -1 );

+mov( U2, OneUlp );

+mov( OneUlp, X );

+for( Indx = 1; Indx <= 3; ++Indx)

+	{

+/*Y = SQRT((X + U1 + X) + F9);*/

+	add( X, U1, Y );

+	add( Y, X, Y );

+	add( Y, F9, Y );

+	mov( Y, Sqarg );

+	SQRT( Sqarg, Y );

+/*Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;*/

+	sub( U2, One, t );

+	add( t, X, t );

+	sub( U2, Y, Y );

+	sub( t, Y, Y );

+	div( OneUlp, Y, Y );

+/*Z = ((U1 - X) + F9) * Half * X * X / OneUlp;*/

+	sub( X, U1, t );

+	add( t, F9, t );

+	mul( t, Half, t );

+	mul( t, X, t );

+	mul( t, X, t );

+	div( OneUlp, t, Z );

+	add( Y, Half, SqEr );

+	add( SqEr, Z, SqEr );

+	Showsq( -1 );

+	sub( Half, Y, SqEr );

+	add( SqEr, Z, SqEr );

+	Showsq( 1 );

+	if(((Indx == 1) || (Indx == 3))) 

+		{

+/*X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));*/

+		mov( OneUlp, Sqarg );

+		SQRT( Sqarg, t );

+		mul( Nine, t, t );

+		div( t, Eight, t );

+		FLOOR( t, t );

+		Sign( X, t2 );

+		mul( t2, t, t );

+		mul( OneUlp, t, X );

+		}

+	else

+		{

+		mov( U1, OneUlp );

+		mov( OneUlp, X );

+		neg( X );

+		}

+	}

+/*=============================================*/

+Milestone = 85;

+/*=============================================*/

+SqRWrng = False;

+Anomaly = False;

+if( cmp(Radix,One) != 0 )

+	{

+	printf("Testing whether sqrt is rounded or chopped.\n");

+/*D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));*/

+	FLOOR( Precision, t2 );

+	add( One, Precision, t );

+	sub( t2, t, t );

+	POW( Radix, t, D );

+	add( Half, D, D );

+	FLOOR( D, D );

+/* ... == Radix^(1 + fract) if (Precision == Integer + fract. */

+	div( Radix, D, X );

+	div( A1, D, Y );

+	FLOOR( X, t );

+	FLOOR( Y, t2 );

+	if( (cmp(X,t) != 0) || (cmp(Y,t2) != 0) )

+		{

+		Anomaly = True;

+		printf( "Anomaly 1\n" );

+		}

+	else

+		{

+		mov( Zero, X );

+		mov( X, Z2 );

+		mov( One, Y );

+		mov( Y, Y2 );

+		sub( One, Radix, Z1 );

+		mul( Four, D, FourD );

+		do

+			{

+			if( cmp(Y2,Z2) >0 )

+				{

+				mov( Radix, Q );

+				mov( Y, YY1 );

+				do

+					{

+/*X1 = FABS(Q + FLOOR(Half - Q / YY1) * YY1);*/

+					div( YY1, Q, t );

+					sub( t, Half, t );

+					FLOOR( t, t );

+					mul( t, YY1, t );

+					add( Q, t, X1 );

+					FABS( X1 );

+					mov( YY1, Q );

+					mov( X1, YY1 );

+					}

+				while( ! (cmp(X1,Zero) <= 0) );

+				if( cmp(Q,One) <= 0 )

+					{

+					mov( Y2, Z2 );

+					mov( Y, Z );

+					}

+				}

+			add( Y, Two, Y );

+			add( X, Eight, X );

+			add( Y2, X, Y2 );

+			if( cmp(Y2,FourD) >= 0 )

+				sub( FourD, Y2, Y2 );

+			}

+		while( ! (cmp(Y,D) >= 0) );

+		sub( Z2, FourD, X8 );

+		mul( Z, Z, Q );

+		add( X8, Q, Q );

+		div( FourD, Q, Q );

+		div( Eight, X8, X8 );

+		FLOOR( Q, t );

+		if( cmp(Q,t) != 0 )

+			{

+			Anomaly = True;

+			printf( "Anomaly 2\n" );

+			}

+		else

+			{

+			Break = False;

+			do

+				{

+				mul( Z1, Z, X );

+/*X = X - FLOOR(X / Radix) * Radix;*/

+				div( Radix, X, t );

+				FLOOR( t, t );

+				mul( t, Radix, t );

+				sub( t, X, X );

+				if( cmp(X,One) == 0 ) 

+					Break = True;

+				else

+					sub( One, Z1, Z1 );

+				}

+			while( ! (Break || (cmp(Z1,Zero) <= 0)) );

+			if( (cmp(Z1,Zero) <= 0) && (! Break))

+				{

+				printf( "Anomaly 3\n" );

+				Anomaly = True;

+				}

+			else

+				{

+				if( cmp(Z1,RadixD2) > 0)

+					sub( Radix, Z1, Z1 );

+				do

+					{

+					NewD();

+					mul( U2, D, t );

+					}

+				while( ! (cmp(t,F9) >= 0) );

+				mul( D, Radix, t );

+				sub( D, t, t );

+				sub( D, W, t2 );

+				if (cmp(t,t2) != 0 )

+					{

+					printf( "Anomaly 4\n" );

+					Anomaly = True;

+					}

+				else

+					{

+					mov( D, Z2 );

+					I = 0;

+					add( One, Z, t );

+					mul( t, Half, t );

+					add( D, t, Y );

+					add( D, Z, X );

+					add( X, Q, X );

+					SR3750();

+					sub( Z, One, t );

+					mul( t, Half, t );

+					add( D, t, Y );

+					add( Y, D, Y );

+					sub( Z, D, X );

+					add( X, D, X );

+					add( X, Q, t );

+					add( t, X, X );

+					SR3750();

+					NewD();

+					sub( Z2, D, t );

+					sub( Z2, W, t2 );

+					if(cmp(t,t2) != 0 )

+						{

+						printf( "Anomaly 5\n" );

+						Anomaly = True;

+						}

+					else

+						{

+/*Y = (D - Z2) + (Z2 + (One - Z) * Half);*/

+						sub( Z, One, t );

+						mul( t, Half, t );

+						add( Z2, t, t );

+						sub( Z2, D, Y );

+						add( Y, t, Y );

+/*X = (D - Z2) + (Z2 - Z + Q);*/

+						sub( Z, Z2, t );

+						add( t, Q, t );

+						sub( Z2, D, X );

+						add( X, t, X );

+						SR3750();

+						add( One, Z, Y );

+						mul( Y, Half, Y );

+						mov( Q, X );

+						SR3750();

+						if(I == 0)

+							{

+							printf( "Anomaly 6\n" );

+							Anomaly = True;

+							}

+						}

+					}

+				}

+			}

+		}

+	if ((I == 0) || Anomaly)

+		{

+		ErrCnt[Failure] += 1;

+		printf( "Anomalous arithmetic with Integer < \n");

+		printf("Radix^Precision = " );

+		show( W );

+		printf(" fails test whether sqrt rounds or chops.\n");

+		SqRWrng = True;

+		}

+	}

+if(! Anomaly)

+	{

+	if(! ((cmp(MinSqEr,Zero) < 0) || (cmp(MaxSqEr,Zero) > 0))) {

+	RSqrt = Rounded;

+	printf("Square root appears to be correctly rounded.\n");

+	}

+	else

+		{

+		k = 0;

+		add( MaxSqEr, U2, t );

+		sub( Half, U2, t2 );

+		if( cmp(t,t2) > 0 )

+			k = 1;

+		if( cmp( MinSqEr, Half ) > 0 )

+			k = 1;

+		add( MinSqEr, Radix, t );

+		if( cmp( t, Half ) < 0 )

+			k = 1;

+		if( k == 1 )

+			SqRWrng = True;

+		else

+			{

+			RSqrt = Chopped;

+			printf("Square root appears to be chopped.\n");

+			}

+		}

+	}

+if( SqRWrng )

+	{

+	printf("Square root is neither chopped nor correctly rounded.\n");

+	printf("Observed errors run from " );

+	sub( Half, MinSqEr, t );

+	show( t );

+	printf("\tto " );

+	add( Half, MaxSqEr, t );

+	show( t );

+	printf( "ulps.\n" );

+	sub( MinSqEr, MaxSqEr, t );

+	mul( Radix, Radix, t2 );

+	if( cmp( t, t2 ) >= 0 )

+		{

+		ErrCnt[Serious] += 1;

+		printf( "sqrt gets too many last digits wrong\n");

+		}

+	}

+}

+

+Showsq( arg )

+int arg;

+{

+

+k = 0;

+if( arg <= 0 )

+	{

+	if( cmp(SqEr,MinSqEr) < 0 )

+		{

+		k = 1;

+		mov( SqEr, MinSqEr );

+		}

+	}

+if( arg >= 0 )

+	{

+	if( cmp(SqEr,MaxSqEr) > 0 )

+		{

+		k = 2;

+		mov( SqEr, MaxSqEr );

+		}

+	}

+#if DEBUG

+if( k != 0 )

+	{

+	printf( "Square root of " );

+	show( arg );

+	printf( "\tis in error by " );

+	show( SqEr );

+	}

+#endif

+}

+

+

+pow1test()

+{

+

+/*=============================================*/

+Milestone = 90;

+/*=============================================*/

+Pause();

+printf("Testing powers Z^i for small Integers Z and i.\n");

+N = 0;

+/* ... test powers of zero. */

+I = 0;

+mov( Zero, Z );

+neg(Z);

+M = 3;

+Break = False;

+do

+	{

+	mov( One, X );

+	SR3980();

+	if(I <= 10)

+		{

+		I = 1023;

+		SR3980();

+		}

+	if( cmp(Z,MinusOne) == 0 )

+		Break = True;

+	else

+		{

+		mov( MinusOne, Z );

+		PrintIfNPositive();

+		N = 0;

+/* .. if(-1)^N is invalid, replace MinusOne by One. */

+		I = - 4;

+		}

+	}

+while( ! Break );

+PrintIfNPositive();

+N1 = N;

+N = 0;

+mov( A1, Z );

+/*M = FLOOR(Two * LOG(W) / LOG(A1));*/

+LOG( W, t );

+mul( Two, t, t );

+FLOOR( t, t );

+LOG( A1, t2 );

+div( t2, t, t );

+FTOL( t, &lngint, t2 );

+M = lngint;

+Break = False;

+do

+	{

+	mov( Z, X );

+	I = 1;

+	SR3980();

+	if( cmp(Z,AInvrse) == 0 )

+		Break = True;

+	else

+		 mov( AInvrse, Z );

+	}

+while( ! (Break) );

+/*=============================================*/

+Milestone = 100;

+/*=============================================*/

+/*  Powers of Radix have been tested, */

+/*         next try a few primes     */

+

+M = NoTrials;

+

+mov( Three, Z );

+do

+	{

+	mov( Z, X );

+	I = 1;

+	SR3980();

+	do

+		{

+		add( Z, Two, Z );

+		div( Three, Z, t );

+		FLOOR( t, t );

+		mul( Three, t, t );

+		}

+	while( cmp(t,Z) == 0 );

+	mul( Eight, Three, t );

+	}

+while( cmp(Z,t) < 0 );

+

+if(N > 0)

+	{

+	printf("Errors like this may invalidate financial calculations\n");

+	printf("\tinvolving interest rates.\n");

+	}

+PrintIfNPositive();

+N += N1;

+if(N == 0)

+	printf("... no discrepancies found.\n");

+if(N > 0)

+	Pause();

+else printf("\n");

+}

+

+

+

+pow2test()

+{

+printf("\n");

+/* ...calculate Exp2 == exp(2) == 7.38905 60989 30650 22723 04275-... */

+mov( Zero, X );

+mov( Two, t2 ); /*I = 2;*/

+

+mul( Two, Three, Y );

+mov( Zero, Q );

+N = 0;

+do

+	{

+	mov( X, Z );

+	add( t2, One, t2 ); /*I = I + 1;*/

+	add( t2, t2, t );

+	div( t, Y, Y ); /*Y = Y / (I + I);*/

+	add( Y, Q, R );

+	add( Z, R, X );

+	sub( X, Z, Q );

+	add( Q, R, Q );

+	}

+while( cmp(X,Z) > 0 );

+

+/*Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);*/

+div( Eight, One, t );

+add( OneAndHalf, t, Z );

+mul( OneAndHalf, ThirtyTwo, t );

+div( t, X, t );

+add( Z, t, Z );

+mul( Z, Z, X );

+mul( X, X, Exp2 );

+mov( F9, X );

+sub( U1, X, Y );

+printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = " );

+show( Exp2 );

+printf( "\tas X -> 1.\n" );

+for(I = 1;;)

+	{

+	sub( BInvrse, X, Z );

+/*Z = (X + One) / (Z - (One - BInvrse));*/

+	add( X, One, t2 );

+	sub( BInvrse, One, t );

+	sub( t, Z, t );

+	div( t, t2, Z );

+	POW( X, Z, Sqarg );

+	sub( Exp2, Sqarg, Q );

+	mov( Q, t );

+	FABS( t );

+	mul( TwoForty, U2, t2 );

+	if( cmp( t, t2 ) > 0 )

+		{

+		N = 1;

+		sub( BInvrse, X, V9 );

+		sub( BInvrse, One, t );

+		sub( t, V9, V9 );

+		ErrCnt[Defect] += 1;

+		printf( "Calculated " );

+		show( Sqarg );

+		printf(" for \t(1 + " );

+		show( V9 );

+		printf( "\tto the power " );

+		show( Z );

+		printf("\tdiffers from correct value by " );

+		show( Q );

+		printf("\tThis much error may spoil financial\n");

+		printf("\tcalculations involving tiny interest rates.\n");

+		break;

+		}

+	else

+		{

+		sub( X, Y, Z );

+		mul( Z, Two, Z );

+		add( Z, Y, Z );

+		mov( Y, X );

+		mov( Z, Y );

+		sub( F9, X, Z );

+		mul( Z, Z, Z );

+		add( Z, One, Z );

+		if( (cmp(Z,One) > 0) && (I < NoTrials) )

+			I++;

+		else

+			{

+			if( cmp(X,One) > 0 )

+				{

+				if(N == 0)

+					printf("Accuracy seems adequate.\n");

+				break;

+				}

+			else

+				{

+				add( One, U2, X );

+				add( U2, U2, Y );

+				add( X, Y, Y );

+				I = 1;

+				}

+			}

+		}

+	}

+/*=============================================*/

+Milestone = 150;

+/*=============================================*/

+printf("Testing powers Z^Q at four nearly extreme values.\n");

+N = 0;

+mov( A1, Z );

+/*Q = FLOOR(Half - LOG(C) / LOG(A1));*/

+LOG( C, t );

+LOG( A1, t2 );

+div( t2, t, t );

+sub( t, Half, t );

+FLOOR( t, Q );

+Break = False;

+do

+	{

+	mov( CInvrse, X );

+	POW( Z, Q, Y );

+	IsYeqX();

+	neg(Q);

+	mov( C, X );

+	POW( Z, Q, Y );

+	IsYeqX();

+	if( cmp(Z,One) < 0 )

+		Break = True;

+	else

+		mov( AInvrse, Z );

+	}

+while( ! (Break));

+PrintIfNPositive();

+if(N == 0)

+	printf(" ... no discrepancies found.\n");

+printf("\n");

+}