run dos2unix on these files

1 files changed, 3550 insertions, 3550 deletions

diff --git a/test/math/eparanoi.c b/test/math/eparanoi.c
index 0e479f9f5..84cab73f8 100644
--- a/test/math/eparanoi.c
+++ b/test/math/eparanoi.c
@@ -1,3550 +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");

-}

+/* 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");
+}