diff options
Diffstat (limited to 'test/math/eparanoi.c')
| -rw-r--r-- | test/math/eparanoi.c | 3550 |
1 files changed, 3550 insertions, 0 deletions
diff --git a/test/math/eparanoi.c b/test/math/eparanoi.c new file mode 100644 index 000000000..0e479f9f5 --- /dev/null +++ b/test/math/eparanoi.c @@ -0,0 +1,3550 @@ +/* paranoia.c arithmetic tester
+ *
+ * This is an implementation of the PARANOIA program. It substitutes
+ * subroutine calls for ALL floating point arithmetic operations.
+ * This permits you to substitute your own experimental versions of
+ * arithmetic routines. It also defeats compiler optimizations,
+ * so for native arithmetic you can be pretty sure you are testing
+ * the arithmetic and not the compiler.
+ *
+ * This version of PARANOIA omits the display of division by zero.
+ * It also omits the test for extra precise subexpressions, since
+ * they cannot occur in this context. Otherwise it includes all the
+ * tests of the 27 Jan 86 distribution, plus a few additional tests.
+ * Commentary has been reduced to a minimum in order to make the program
+ * smaller.
+ *
+ * The original PARANOIA program, written by W. Kahan, C version
+ * by Thos Sumner and David Gay, can be downloaded free from the
+ * Internet NETLIB. An MSDOS disk can be obtained for $15 from:
+ * Richard Karpinski
+ * 6521 Raymond Street
+ * Oakland, CA 94609
+ *
+ * Steve Moshier, 28 Oct 88
+ * last rev: 23 May 92
+ */
+
+#define DEBUG 0
+
+/* To use the native arithmetic of the computer, define NATIVE
+ * to be 1. To use your own supplied arithmetic routines, NATIVE is 0.
+ */
+#define NATIVE 0
+
+/* gcc real.c interface */
+#define L128DOUBLE 0
+
+#include <stdio.h>
+
+
+
+
+/* Data structure of floating point number. If NATIVE was
+ * selected above, you can define LDOUBLE 1 to test 80-bit long double
+ * precision or define it 0 to test 64-bit double precision.
+*/
+#define LDOUBLE 0
+#if NATIVE
+
+#define NE 1
+#if LDOUBLE
+#define FSIZE long double
+#define FLOAT(x) FSIZE x[NE]
+static FSIZE eone[NE] = {1.0L}; /* The constant 1.0 */
+#define ZSQRT sqrtl
+#define ZLOG logl
+#define ZFLOOR floorl
+#define ZPOW powl
+long double sqrtl(), logl(), floorl(), powl();
+#define FSETUP einit
+#else /* not LDOUBLE */
+#define FSIZE double
+#define FLOAT(x) FSIZE x[NE]
+static FSIZE eone[NE] = {1.0}; /* The constant 1.0 */
+#define ZSQRT sqrt
+#define ZLOG log
+#define ZFLOOR floor
+#define ZPOW pow
+double sqrt(), log(), floor(), pow();
+/* Coprocessor initialization,
+ * defeat underflow trap or what have you.
+ * This is required mainly on i386 and 68K processors.
+ */
+#define FSETUP dprec
+#endif /* double, not LDOUBLE */
+
+#else /* not NATIVE */
+
+/* Setup for extended double type.
+ * Put NE = 10 for real.c operating with TFmode support (16-byte reals)
+ * Put NE = 6 for real.c operating with XFmode support (10- or 12-byte reals)
+ * The value of NE must agree with that in ehead.h, if ieee.c is used.
+ */
+#define NE 6
+#define FSIZE unsigned short
+#define FLOAT(x) unsigned short x[NE]
+extern unsigned short eone[];
+#define FSETUP einit
+
+/* default for FSETUP */
+/*
+einit()
+{}
+*/
+
+error(s)
+char *s;
+{
+printf( "error: %s\n", s );
+}
+
+#endif /* not NATIVE */
+
+
+
+#if L128DOUBLE
+/* real.c interface */
+
+#undef FSETUP
+#define FSETUP efsetup
+
+FLOAT(enone);
+
+#define ONE enone
+
+/* Use emov to convert from widest type to widest type, ... */
+/*
+#define ENTOE emov
+#define ETOEN emov
+*/
+
+/* ... else choose e24toe, e53toe, etc. */
+#define ENTOE e64toe
+#define ETOEN etoe64
+#define NNBITS 64
+
+#define NIBITS ((NE-1)*16)
+extern int rndprc;
+
+efsetup()
+{
+rndprc = NNBITS;
+ETOEN(eone, enone);
+}
+
+add(a,b,c)
+FLOAT(a);
+FLOAT(b);
+FLOAT(c);
+{
+unsigned short aa[10], bb[10], cc[10];
+
+ENTOE(a,aa);
+ENTOE(b,bb);
+eadd(aa,bb,cc);
+ETOEN(cc,c);
+}
+
+sub(a,b,c)
+FLOAT(a);
+FLOAT(b);
+FLOAT(c);
+{
+unsigned short aa[10], bb[10], cc[10];
+
+ENTOE(a,aa);
+ENTOE(b,bb);
+esub(aa,bb,cc);
+ETOEN(cc,c);
+}
+
+mul(a,b,c)
+FLOAT(a);
+FLOAT(b);
+FLOAT(c);
+{
+unsigned short aa[10], bb[10], cc[10];
+
+ENTOE(a,aa);
+ENTOE(b,bb);
+emul(aa,bb,cc);
+ETOEN(cc,c);
+}
+
+div(a,b,c)
+FLOAT(a);
+FLOAT(b);
+FLOAT(c);
+{
+unsigned short aa[10], bb[10], cc[10];
+
+ENTOE(a,aa);
+ENTOE(b,bb);
+ediv(aa,bb,cc);
+ETOEN(cc,c);
+}
+
+int cmp(a,b)
+FLOAT(a);
+FLOAT(b);
+{
+unsigned short aa[10], bb[10];
+int c;
+int ecmp();
+
+ENTOE(a,aa);
+ENTOE(b,bb);
+c = ecmp(aa,bb);
+return(c);
+}
+
+mov(a,b)
+FLOAT(a);
+FLOAT(b);
+{
+int i;
+
+for( i=0; i<NE; i++ )
+ b[i] = a[i];
+}
+
+
+neg(a)
+FLOAT(a);
+{
+unsigned short aa[10];
+
+ENTOE(a,aa);
+eneg(aa);
+ETOEN(aa,a);
+}
+
+clear(a)
+FLOAT(a);
+{
+int i;
+
+for( i=0; i<NE; i++ )
+ a[i] = 0;
+}
+
+FABS(a)
+FLOAT(a);
+{
+unsigned short aa[10];
+
+ENTOE(a,aa);
+eabs(aa);
+ETOEN(aa,a);
+}
+
+FLOOR(a,b)
+FLOAT(a);
+FLOAT(b);
+{
+unsigned short aa[10], bb[10];
+
+ENTOE(a,aa);
+efloor(aa,bb);
+ETOEN(bb,b);
+}
+
+LOG(a,b)
+FLOAT(a);
+FLOAT(b);
+{
+unsigned short aa[10], bb[10];
+int rndsav;
+
+ENTOE(a,aa);
+rndsav = rndprc;
+rndprc = NIBITS;
+elog(aa,bb);
+rndprc = rndsav;
+ETOEN(bb,b);
+}
+
+POW(a,b,c)
+FLOAT(a);
+FLOAT(b);
+FLOAT(c);
+{
+unsigned short aa[10], bb[10], cc[10];
+int rndsav;
+
+ENTOE(a,aa);
+ENTOE(b,bb);
+rndsav = rndprc;
+rndprc = NIBITS;
+epow(aa,bb,cc);
+rndprc = rndsav;
+ETOEN(cc,c);
+}
+
+SQRT(a,b)
+FLOAT(a);
+FLOAT(b);
+{
+unsigned short aa[10], bb[10];
+
+ENTOE(a,aa);
+esqrt(aa,bb);
+ETOEN(bb,b);
+}
+
+FTOL(x,ip,f)
+FLOAT(x);
+long *ip;
+FLOAT(f);
+{
+unsigned short xx[10], ff[10];
+
+ENTOE(x,xx);
+eifrac(xx,ip,ff);
+ETOEN(ff,f);
+}
+
+LTOF(ip,x)
+long *ip;
+FLOAT(x);
+{
+unsigned short xx[10];
+ltoe(ip,xx);
+ETOEN(xx,x);
+}
+
+TOASC(a,b,c)
+FLOAT(a);
+int b;
+char *c;
+{
+unsigned short xx[10];
+
+ENTOE(a,xx);
+etoasc(xx,b,c);
+}
+
+#else /* not L128DOUBLE */
+
+#define ONE eone
+
+/* Note all arguments of operation subroutines are pointers. */
+/* c = b + a */
+#define add(a,b,c) eadd(a,b,c)
+/* c = b - a */
+#define sub(a,b,c) esub(a,b,c)
+/* c = b * a */
+#define mul(a,b,c) emul(a,b,c)
+/* c = b / a */
+#define div(a,b,c) ediv(a,b,c)
+/* 1 if a>b, 0 if a==b, -1 if a<b */
+#define cmp(a,b) ecmp(a,b)
+/* b = a */
+#define mov(a,b) emov(a,b)
+/* a = -a */
+#define neg(a) eneg(a)
+/* a = 0 */
+#define clear(a) eclear(a)
+
+#define FABS(x) eabs(x)
+#define FLOOR(x,y) efloor(x,y)
+#define LOG(x,y) elog(x,y)
+#define POW(x,y,z) epow(x,y,z)
+#define SQRT(x,y) esqrt(x,y)
+
+/* x = &FLOAT input, i = &long integer part, f = &FLOAT fractional part */
+#define FTOL(x,i,f) eifrac(x,i,f)
+
+/* i = &long integer input, x = &FLOAT output */
+#define LTOF(i,x) ltoe(i,x)
+
+/* Convert FLOAT a to decimal ASCII string with b digits */
+#define TOASC(a,b,c) etoasc(a,b,c)
+#endif /* not L128DOUBLE */
+
+
+
+/* The following subroutines are implementations of the above
+ * named functions, using the native or default arithmetic.
+ */
+#if NATIVE
+eadd(a,b,c)
+FSIZE *a, *b, *c;
+{
+*c = *b + *a;
+}
+
+esub(a,b,c)
+FSIZE *a, *b, *c;
+{
+*c = *b - *a;
+}
+
+emul(a,b,c)
+FSIZE *a, *b, *c;
+{
+*c = (*b) * (*a);
+}
+
+ediv(a,b,c)
+FSIZE *a, *b, *c;
+{
+*c = (*b) / (*a);
+}
+
+
+/* Important note: comparison can be done by subracting
+ * or by a compare instruction that may or may not be
+ * equivalent to subtracting.
+ */
+ecmp(a,b)
+FSIZE *a, *b;
+{
+if( (*a) > (*b) )
+ return( 1 );
+if( (*a) < (*b) )
+ return( -1 );
+if( (*a) != (*b) )
+ goto cmpf;
+if( (*a) == (*b) )
+ return( 0 );
+cmpf:
+printf( "Compare fails\n" );
+return(0);
+}
+
+
+emov( a, b )
+FSIZE *a, *b;
+{
+*b = *a;
+}
+
+eneg( a )
+FSIZE *a;
+{
+*a = -(*a);
+}
+
+eclear(a)
+FSIZE *a;
+{
+*a = 0.0;
+}
+
+eabs(x)
+FSIZE *x;
+{
+if( (*x) < 0.0 )
+ *x = -(*x);
+}
+
+efloor(x,y)
+FSIZE *x, *y;
+{
+
+*y = (FSIZE )ZFLOOR( *x );
+}
+
+elog(x,y)
+FSIZE *x, *y;
+{
+
+*y = (FSIZE )ZLOG( *x );
+}
+
+epow(x,y,z)
+FSIZE *x, *y, *z;
+{
+
+*z = (FSIZE )ZPOW( *x, *y );
+}
+
+esqrt(x,y)
+FSIZE *x, *y;
+{
+
+*y = (FSIZE )ZSQRT( *x );
+}
+
+
+eifrac(x,i,f)
+FSIZE *x;
+long *i;
+FSIZE *f;
+{
+FSIZE y;
+
+y = (FSIZE )ZFLOOR( *x );
+if( y < 0.0 )
+ {
+ *f = y - *x;
+ *i = -y;
+ }
+else
+ {
+ *f = *x - y;
+ *i = y;
+ }
+}
+
+
+ltoe(i,x)
+long *i;
+FSIZE *x;
+{
+*x = *i;
+}
+
+
+etoasc(a,str,n)
+FSIZE *a;
+char *str;
+int n;
+{
+double x;
+
+x = (double )(*a);
+sprintf( str, " %.17e ", x );
+}
+
+/* default for FSETUP */
+einit()
+{}
+
+#endif /* NATIVE */
+
+
+
+
+FLOAT(Radix);
+FLOAT(BInvrse);
+FLOAT(RadixD2);
+FLOAT(BMinusU2);
+/*Small floating point constants.*/
+FLOAT(Zero);
+FLOAT(Half);
+FLOAT(One);
+FLOAT(Two);
+FLOAT(Three);
+FLOAT(Four);
+FLOAT(Five);
+FLOAT(Six);
+FLOAT(Eight);
+FLOAT(Nine);
+FLOAT(Ten);
+FLOAT(TwentySeven);
+FLOAT(ThirtyTwo);
+FLOAT(TwoForty);
+FLOAT(MinusOne );
+FLOAT(OneAndHalf);
+
+/*Integer constants*/
+int NoTrials = 20; /*Number of tests for commutativity. */
+#define False 0
+#define True 1
+
+/* Definitions for declared types
+ Guard == (Yes, No);
+ Rounding == (Chopped, Rounded, Other);
+ Message == packed array [1..40] of char;
+ Class == (Flaw, Defect, Serious, Failure);
+ */
+#define Yes 1
+#define No 0
+#define Chopped 2
+#define Rounded 1
+#define Other 0
+#define Flaw 3
+#define Defect 2
+#define Serious 1
+#define Failure 0
+
+typedef int Guard, Rounding, Class;
+typedef char Message;
+
+/* Declarations of Variables */
+FLOAT(AInvrse);
+FLOAT(A1);
+FLOAT(C);
+FLOAT(CInvrse);
+FLOAT(D);
+FLOAT(FourD);
+FLOAT(E0);
+FLOAT(E1);
+FLOAT(Exp2);
+FLOAT(E3);
+FLOAT(MinSqEr);
+FLOAT(SqEr);
+FLOAT(MaxSqEr);
+FLOAT(E9);
+FLOAT(Third);
+FLOAT(F6);
+FLOAT(F9);
+FLOAT(H);
+FLOAT(HInvrse);
+FLOAT(StickyBit);
+FLOAT(J);
+FLOAT(MyZero);
+FLOAT(Precision);
+FLOAT(Q);
+FLOAT(Q9);
+FLOAT(R);
+FLOAT(Random9);
+FLOAT(T);
+FLOAT(Underflow);
+FLOAT(S);
+FLOAT(OneUlp);
+FLOAT(UfThold);
+FLOAT(U1);
+FLOAT(U2);
+FLOAT(V);
+FLOAT(V0);
+FLOAT(V9);
+FLOAT(W);
+FLOAT(X);
+FLOAT(X1);
+FLOAT(X2);
+FLOAT(X8);
+FLOAT(Random1);
+FLOAT(Y);
+FLOAT(YY1);
+FLOAT(Y2);
+FLOAT(Random2);
+FLOAT(Z);
+FLOAT(PseudoZero);
+FLOAT(Z1);
+FLOAT(Z2);
+FLOAT(Z9);
+static FLOAT(t);
+FLOAT(t2);
+FLOAT(Sqarg);
+int ErrCnt[4];
+int fpecount;
+int Milestone;
+int PageNo;
+int I, M, N, N1, stkflg;
+Guard GMult, GDiv, GAddSub;
+Rounding RMult, RDiv, RAddSub, RSqrt;
+int Break, Done, NotMonot, Monot, Anomaly, IEEE;
+int SqRWrng, UfNGrad;
+int k, k2;
+int Indx;
+char ch[8];
+
+long lngint, lng2; /* intermediate for conversion between int and FLOAT */
+
+/* Computed constants. */
+/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
+/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
+
+
+show( x )
+short x[];
+{
+int i;
+char s[80];
+
+/* Number of 16-bit groups to display */
+#if NATIVE
+#if LDOUBLE
+#define NPRT (sizeof( long double )/2)
+#else
+#define NPRT (sizeof( double )/2)
+#endif
+#else
+#define NPRT NE
+#endif
+
+TOASC( x, s, 70 );
+printf( "%s\n", s );
+for( i=0; i<NPRT; i++ )
+ printf( "%04x ", x[i] & 0xffff );
+printf( "\n" );
+}
+
+/* define NOSIGNAL */
+#ifndef NOSIGNAL
+#include <signal.h>
+#endif
+#include <setjmp.h>
+jmp_buf ovfl_buf;
+/*typedef int (*Sig_type)();*/
+typedef void (*Sig_type)();
+Sig_type sigsave;
+
+/* Floating point exception receiver */
+void sigfpe()
+{
+fpecount++;
+printf( "\n* * * FLOATING-POINT ERROR * * *\n" );
+/* reinitialize the floating point unit */
+FSETUP();
+fflush(stdout);
+if( sigsave )
+ {
+#ifndef NOSIGNAL
+ signal( SIGFPE, sigsave );
+#endif
+ sigsave = 0;
+ longjmp( ovfl_buf, 1 );
+ }
+abort();
+}
+
+
+main()
+{
+
+/* Do coprocessor or other initializations */
+FSETUP();
+
+printf(
+ "This version of paranoia omits test for extra precise subexpressions\n" );
+printf( "and includes a few additional tests.\n" );
+
+clear(Zero);
+printf( "0 = " );
+show( Zero );
+mov( ONE, One);
+printf( "1 = " );
+show( One );
+add( One, One, Two );
+printf( "1+1 = " );
+show( Two );
+add( Two, One, Three );
+add( Three, One, Four );
+add( Four, One, Five );
+add( Five, One, Six );
+add( Four, Four, Eight );
+mul( Three, Three, Nine );
+add( Nine, One, Ten );
+mul( Nine, Three, TwentySeven );
+mul( Four, Eight, ThirtyTwo );
+mul( Four, Five, t );
+mul( t, Three, t );
+mul( t, Four, TwoForty );
+mov( One, MinusOne );
+neg( MinusOne );
+div( Two, One, Half );
+add( One, Half, OneAndHalf );
+ErrCnt[Failure] = 0;
+ErrCnt[Serious] = 0;
+ErrCnt[Defect] = 0;
+ErrCnt[Flaw] = 0;
+PageNo = 1;
+#ifndef NOSIGNAL
+signal( SIGFPE, sigfpe );
+#endif
+printf("Program is now RUNNING tests on small integers:\n");
+
+add( Zero, Zero, t );
+if( cmp( t, Zero ) != 0)
+ {
+ printf( "0+0 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+sub( One, One, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "1-1 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+if( cmp( One, Zero ) <= 0 )
+ {
+ printf( "1 <= 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( One, One, t );
+if( cmp( t, Two ) != 0 )
+ {
+ printf( "1+1 != 2\n" );
+ ErrCnt[Failure] += 1;
+ }
+mov( Zero, Z );
+neg( Z );
+FLOOR( Z, t );
+if( cmp(t,Zero) != 0 )
+ {
+ ErrCnt[Serious] += 1;
+ printf( "FLOOR(-0) should equal 0, is = " );
+ show( t );
+ }
+if( cmp(Z, Zero) != 0)
+ {
+ ErrCnt[Failure] += 1;
+ printf("Comparison alleges that -0.0 is Non-zero!\n");
+ }
+else
+ {
+ div( TwoForty, One, U1 ); /* U1 = 0.001 */
+ mov( One, Radix );
+ TstPtUf();
+ }
+add( Two, One, t );
+if( cmp( t, Three ) != 0 )
+ {
+ printf( "2+1 != 3\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( Three, One, t );
+if( cmp( t, Four ) != 0 )
+ {
+ printf( "3+1 != 4\n" );
+ ErrCnt[Failure] += 1;
+ }
+mov( Two, t );
+neg( t );
+mul( Two, t, t );
+add( Four, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "4+2*(-2) != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+sub( Three, Four, t );
+sub( One, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "4-3-1 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+ sub( One, Zero, t );
+if( cmp( t, MinusOne ) != 0 )
+ {
+ printf( "-1 != 0-1\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( One, MinusOne, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "1+(-1) != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+mov( One, t );
+FABS( t );
+add( MinusOne, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "-1+abs(1) != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+mul( MinusOne, MinusOne, t );
+add( MinusOne, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "-1+(-1)*(-1) != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( Half, MinusOne, t );
+add( Half, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "1/2 + (-1) + 1/2 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+Milestone = 10;
+mul( Three, Three, t );
+if( cmp( t, Nine ) != 0 )
+ {
+ printf( "3*3 != 9\n" );
+ ErrCnt[Failure] += 1;
+ }
+mul( Nine, Three, t );
+if( cmp( t, TwentySeven ) != 0 )
+ {
+ printf( "3*9 != 27\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( Four, Four, t );
+if( cmp( t, Eight ) != 0 )
+ {
+ printf( "4+4 != 8\n" );
+ ErrCnt[Failure] += 1;
+ }
+mul( Eight, Four, t );
+if( cmp( t, ThirtyTwo ) != 0 )
+ {
+ printf( "8*4 != 32\n" );
+ ErrCnt[Failure] += 1;
+ }
+sub( TwentySeven, ThirtyTwo, t );
+sub( Four, t, t );
+sub( One, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "32-27-4-1 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( Four, One, t );
+if( cmp( t, Five ) != 0 )
+ {
+ printf( "4+1 != 5\n" );
+ ErrCnt[Failure] += 1;
+ }
+mul( Four, Five, t );
+mul( Three, t, t );
+mul( Four, t, t );
+if( cmp( t, TwoForty ) != 0 )
+ {
+ printf( "4*5*3*4 != 240\n" );
+ ErrCnt[Failure] += 1;
+ }
+div( Three, TwoForty, t );
+mul( Four, Four, t2 );
+mul( Five, t2, t2 );
+sub( t2, t2, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "240/3 - 4*4*5 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+div( Four, TwoForty, t );
+mul( Five, Three, t2 );
+mul( Four, t2, t2 );
+sub( t2, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "240/4 - 5*3*4 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+div( Five, TwoForty, t );
+mul( Four, Three, t2 );
+mul( Four, t2, t2 );
+sub( t2, t, t );
+if( cmp( t, Zero ) != 0 )
+ {
+ printf( "240/5 - 4*3*4 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+if(ErrCnt[Failure] == 0)
+ {
+printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n\n");
+ }
+printf("Searching for Radix and Precision.\n");
+mov( One, W );
+do
+ {
+ add( W, W, W );
+ add( W, One, Y );
+ sub( W, Y, Z );
+ sub( One, Z, Y );
+ mov( Y, t );
+ FABS(t);
+ add( MinusOne, t, t );
+ k = cmp( t, Zero );
+ }
+while( k < 0 );
+/*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
+mov( Zero, Precision );
+mov( One, Y );
+do
+ {
+ add( W, Y, Radix );
+ add( Y, Y, Y );
+ sub( W, Radix, Radix );
+ k = cmp( Radix, Zero );
+ }
+while( k == 0);
+
+if( cmp(Radix, Two) < 0 )
+ mov( One, Radix );
+printf("Radix = " );
+show( Radix );
+if( cmp(Radix, One) != 0)
+ {
+ mov( One, W );
+ do
+ {
+ add( One, Precision, Precision );
+ mul( W, Radix, W );
+ add( W, One, Y );
+ sub( W, Y, t );
+ k = cmp( t, One );
+ }
+ while( k == 0 );
+ }
+/* now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 */
+div( W, One, U1 );
+mul( Radix, U1, U2 );
+printf( "Closest relative separation found is U 1 = " );
+show( U1 );
+printf( "Recalculating radix and precision." );
+
+/*save old values*/
+mov( Radix, E0 );
+mov( U1, E1 );
+mov( U2, E9 );
+mov( Precision, E3 );
+
+div( Three, Four, X );
+sub( One, X, Third );
+sub( Third, Half, F6 );
+add( F6, F6, X );
+sub( Third, X, X );
+FABS( X );
+if( cmp(X, U2) < 0 )
+ mov( U2, X );
+
+/*... now X = (unknown no.) ulps of 1+...*/
+do
+ {
+ mov( X, U2 );
+/* Y = Half * U2 + ThirtyTwo * U2 * U2; */
+ mul( ThirtyTwo, U2, t );
+ mul( t, U2, t );
+ mul( Half, U2, Y );
+ add( t, Y, Y );
+ add( One, Y, Y );
+ sub( One, Y, X );
+ k = cmp( U2, X );
+ k2 = cmp( X, Zero );
+ }
+while ( ! ((k <= 0) || (k2 <= 0)));
+
+/*... now U2 == 1 ulp of 1 + ... */
+div( Three, Two, X );
+sub( Half, X, F6 );
+add( F6, F6, Third );
+sub( Half, Third, X );
+add( F6, X, X );
+FABS( X );
+if( cmp(X, U1) < 0 )
+ mov( U1, X );
+
+/*... now X == (unknown no.) ulps of 1 -... */
+do
+ {
+ mov( X, U1 );
+ /* Y = Half * U1 + ThirtyTwo * U1 * U1;*/
+ mul( ThirtyTwo, U1, t );
+ mul( U1, t, t );
+ mul( Half, U1, Y );
+ add( t, Y, Y );
+ sub( Y, Half, Y );
+ add( Half, Y, X );
+ sub( X, Half, Y );
+ add( Half, Y, X );
+ k = cmp( U1, X );
+ k2 = cmp( X, Zero );
+ } while ( ! ((k <= 0) || (k2 <= 0)));
+/*... now U1 == 1 ulp of 1 - ... */
+if( cmp( U1, E1 ) == 0 )
+ printf("confirms closest relative separation U1 .\n");
+else
+ {
+ printf("gets better closest relative separation U1 = " );
+ show( U1 );
+ }
+div( U1, One, W );
+sub( U1, Half, F9 );
+add( F9, Half, F9 );
+div( U1, U2, t );
+div( TwoForty, One, t2 );
+add( t2, t, t );
+FLOOR( t, Radix );
+if( cmp(Radix, E0) == 0 )
+ printf("Radix confirmed.\n");
+else
+ {
+ printf("MYSTERY: recalculated Radix = " );
+ show( Radix );
+ mov( E0, Radix );
+ }
+add( Eight, Eight, t );
+if( cmp( Radix, t ) > 0 )
+ {
+ printf( "Radix is too big: roundoff problems\n" );
+ ErrCnt[Defect] += 1;
+ }
+k = 1;
+if( cmp( Radix, Two ) == 0 )
+ k = 0;
+if( cmp( Radix, Ten ) == 0 )
+ k = 0;
+if( cmp( Radix, One ) == 0 )
+ k = 0;
+if( k != 0 )
+ {
+ printf( "Radix is not as good as 2 or 10\n" );
+ ErrCnt[Flaw] += 1;
+ }
+/*=============================================*/
+Milestone = 20;
+/*=============================================*/
+sub( Half, F9, t );
+if( cmp( t, Half ) >= 0 )
+ {
+ printf( "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?\n" );
+ ErrCnt[Failure] += 1;
+ }
+mov( F9, X );
+I = 1;
+sub( Half, X, Y );
+sub( Half, Y, Z );
+if( (cmp( X, One ) == 0) && (cmp( Z, Zero) != 0) )
+ {
+ printf( "Comparison is fuzzy ,X=1 but X-1/2-1/2 != 0\n" );
+ ErrCnt[Failure] += 1;
+ }
+add( One, U2, X );
+I = 0;
+/*=============================================*/
+Milestone = 25;
+/*=============================================*/
+/*... BMinusU2 = nextafter(Radix, 0) */
+
+sub( One, Radix, BMinusU2 );
+sub( U2, BMinusU2, t );
+add( One, t, BMinusU2 );
+/* Purify Integers */
+if( cmp(Radix,One) != 0 )
+ {
+/*X = - TwoForty * LOG(U1) / LOG(Radix);*/
+ LOG( U1, X );
+ LOG( Radix, t );
+ div( t, X, X );
+ mul( TwoForty, X, X );
+ neg( X );
+
+ add( Half, X, Y );
+ FLOOR( Y, Y );
+ sub( Y, X, t );
+ FABS( t );
+ mul( Four, t, t );
+ if( cmp( t, One ) < 0 )
+ mov( Y, X );
+ div( TwoForty, X, Precision );
+ add( Half, Precision, Y );
+ FLOOR( Y, Y );
+ sub( Y, Precision, t );
+ FABS( t );
+ mul( TwoForty, t, t );
+ if( cmp( t, Half ) < 0 )
+ mov( Y, Precision );
+ }
+FLOOR( Precision, t );
+if( (cmp( Precision, t ) != 0) || (cmp( Radix, One ) == 0) )
+ {
+ printf("Precision cannot be characterized by an Integer number\n");
+ printf("of significant digits but, by itself, this is a minor flaw.\n");
+ }
+if( cmp(Radix, One) == 0 )
+ printf("logarithmic encoding has precision characterized solely by U1.\n");
+else
+ {
+ printf("The number of significant digits of the Radix is " );
+ show( Precision );
+ }
+mul( U2, Nine, t );
+mul( Nine, t, t );
+mul( TwoForty, t, t );
+if( cmp( t, One ) >= 0 )
+ {
+ printf( "Precision worse than 5 decimal figures\n" );
+ ErrCnt[Serious] += 1;
+ }
+/*=============================================*/
+Milestone = 30;
+/*=============================================*/
+/* Test for extra-precise subepressions has been deleted. */
+Milestone = 35;
+/*=============================================*/
+if( cmp(Radix,Two) >= 0 )
+ {
+ mul( Radix, Radix, t );
+ div( t, W, X );
+ add( X, One, Y );
+ sub( X, Y, Z );
+ add( Z, U2, T );
+ sub( Z, T, X );
+ if( cmp( X, U2 ) != 0 )
+ {
+ printf( "Subtraction is not normalized X=Y,X+Z != Y+Z!\n" );
+ ErrCnt[Failure] += 1;
+ }
+ if( cmp(X,U2) == 0 )
+ printf("Subtraction appears to be normalized, as it should be.");
+ }
+
+printf("\nChecking for guard digit in *, /, and -.\n");
+mul( F9, One, Y );
+mul( One, F9, Z );
+sub( Half, F9, X );
+sub( Half, Y, Y );
+sub( X, Y, Y );
+sub( Half, Z, Z );
+sub( X, Z, Z );
+add( One, U2, X );
+mul( X, Radix, T );
+mul( Radix, X, R );
+sub( Radix, T, X );
+mul( Radix, U2, t );
+sub( t, X, X );
+sub( Radix, R, T );
+mul( Radix, U2, t );
+sub( t, T, T );
+sub( One, Radix, t );
+mul( t, X, X );
+sub( One, Radix, t );
+mul( t, T, T );
+
+k = cmp(X,Zero);
+k |= cmp(Y,Zero);
+k |= cmp(Z,Zero);
+k |= cmp(T,Zero);
+if( k == 0 )
+ GMult = Yes;
+else
+ {
+ GMult = No;
+ ErrCnt[Serious] += 1;
+ printf( "* lacks a Guard Digit, so 1*X != X\n" );
+ }
+mul( Radix, U2, Z );
+add( One, Z, X );
+add( X, Z, Y );
+mul( X, X, t );
+sub( t, Y, Y );
+FABS( Y );
+sub( U2, Y, Y );
+sub( U2, One, X );
+sub( U2, X, Z );
+mul( X, X, t );
+sub( t, Z, Z );
+FABS( Z );
+sub( U1, Z, Z );
+if( (cmp(Y,Zero) > 0) || (cmp(Z,Zero) > 0) )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "* gets too many final digits wrong.\n" );
+ }
+sub( U2, One, Y );
+add( One, U2, X );
+div( Y, One, Z );
+sub( X, Z, Y );
+div( Three, One, X );
+div( Nine, Three, Z );
+sub( Z, X, X );
+div( TwentySeven, Nine, T );
+sub( T, Z, Z );
+k = cmp( X, Zero );
+k |= cmp( Y, Zero );
+k |= cmp( Z, Zero );
+if( k )
+ {
+ ErrCnt[Defect] += 1;
+printf( "Division lacks a Guard Digit, so error can exceed 1 ulp\n" );
+printf( "or 1/3 and 3/9 and 9/27 may disagree\n" );
+ }
+div( One, F9, Y );
+sub( Half, F9, X );
+sub( Half, Y, Y );
+sub( X, Y, Y );
+add( One, U2, X );
+div( One, X, T );
+sub( X, T, X );
+k = cmp( X, Zero );
+k |= cmp( Y, Zero );
+k |= cmp( Z, Zero );
+if( k == 0 )
+ GDiv = Yes;
+else
+ {
+ GDiv = No;
+ ErrCnt[Serious] += 1;
+ printf( "Division lacks a Guard Digit, so X/1 != X\n" );
+ }
+add( One, U2, X );
+div( X, One, X );
+sub( Half, X, Y );
+sub( Half, Y, Y );
+if( cmp(Y,Zero) >= 0 )
+ {
+ ErrCnt[Serious] += 1;
+ printf( "Computed value of 1/1.000..1 >= 1\n" );
+ }
+sub( U2, One, X );
+mul( Radix, U2, Y );
+add( One, Y, Y );
+mul( X, Radix, Z );
+mul( Y, Radix, T );
+div( Radix, Z, R );
+div( Radix, T, StickyBit );
+sub( X, R, X );
+sub( Y, StickyBit, Y );
+k = cmp( X, Zero );
+k |= cmp( Y, Zero );
+if( k )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "* and/or / gets too many last digits wrong\n" );
+ }
+sub( U1, One, Y );
+sub( F9, One, X );
+sub( Y, One, Y );
+sub( U2, Radix, T );
+sub( BMinusU2, Radix, Z );
+sub( T, Radix, T );
+k = cmp( X, U1 );
+k |= cmp( Y, U1 );
+k |= cmp( Z, U2 );
+k |= cmp( T, U2 );
+if( k == 0 )
+ GAddSub = Yes;
+else
+ {
+ GAddSub = No;
+ ErrCnt[Serious] += 1;
+ printf( "- lacks Guard Digit, so cancellation is obscured\n" );
+ }
+sub( One, F9, t );
+if( (cmp(F9,One) != 0) && (cmp(t,Zero) >= 0) )
+ {
+ ErrCnt[Serious] += 1;
+ printf("comparison alleges (1-U1) < 1 although\n");
+ printf(" subtration yields (1-U1) - 1 = 0 , thereby vitiating\n");
+ printf(" such precautions against division by zero as\n");
+ printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
+ }
+if (GMult == Yes && GDiv == Yes && GAddSub == Yes)
+ printf(" *, /, and - appear to have guard digits, as they should.\n");
+/*=============================================*/
+Milestone = 40;
+/*=============================================*/
+printf("Checking rounding on multiply, divide and add/subtract.\n");
+RMult = Other;
+RDiv = Other;
+RAddSub = Other;
+div( Two, Radix, RadixD2 );
+mov( Two, A1 );
+Done = False;
+do
+ {
+ mov( Radix, AInvrse );
+ do
+ {
+ mov( AInvrse, X );
+ div( A1, AInvrse, AInvrse );
+ FLOOR( AInvrse, t );
+ k = cmp( t, AInvrse );
+ }
+ while( ! (k != 0 ) );
+ k = cmp( X, One );
+ k2 = cmp( A1, Three );
+ Done = (k == 0) || (k2 > 0);
+ if(! Done)
+ add( Nine, One, A1 );
+ }
+while( ! (Done));
+if( cmp(X, One) == 0 )
+ mov( Radix, A1 );
+div( A1, One, AInvrse );
+mov( A1, X );
+mov( AInvrse, Y );
+Done = False;
+do
+ {
+ mul( X, Y, Z );
+ sub( Half, Z, Z );
+ if( cmp( Z, Half ) != 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "X * (1/X) differs from 1\n" );
+ }
+ k = cmp( X, Radix );
+ Done = (k == 0);
+ mov( Radix, X );
+ div( X, One, Y );
+ }
+while( ! (Done));
+
+add( One, U2, Y2 );
+sub( U2, One, YY1 );
+sub( U2, OneAndHalf, X );
+add( OneAndHalf, U2, Y );
+sub( U2, X, Z );
+mul( Z, Y2, Z );
+mul( Y, YY1, T );
+sub( X, Z, Z );
+sub( X, T, T );
+mul( X, Y2, X );
+add( Y, U2, Y );
+mul( Y, YY1, Y );
+sub( OneAndHalf, X, X );
+sub( OneAndHalf, Y, Y );
+k = cmp( X, Zero );
+k |= cmp( Y, Zero );
+k |= cmp( Z, Zero );
+if( cmp( T, Zero ) > 0 )
+ k = 1;
+if( k == 0 )
+ {
+ add( OneAndHalf, U2, X );
+ mul( X, Y2, X );
+ sub( U2, OneAndHalf, Y );
+ sub( U2, Y, Y );
+ add( OneAndHalf, U2, Z );
+ add( U2, Z, Z );
+ sub( U2, OneAndHalf, T );
+ mul( T, YY1, T );
+ add( Z, U2, t );
+ sub( t, X, X );
+ mul( Y, YY1, StickyBit );
+ mul( Z, Y2, S );
+ sub( Y, T, T );
+ sub( Y, U2, Y );
+ add( StickyBit, Y, Y );
+/* Z = S - (Z + U2 + U2); */
+ add( Z, U2, t );
+ add( t, U2, t );
+ sub( t, S, Z );
+ add( Y2, U2, t );
+ mul( t, YY1, StickyBit );
+ mul( Y2, YY1, YY1 );
+ sub( Y2, StickyBit, StickyBit );
+ sub( Half, YY1, YY1 );
+ k = cmp( X, Zero );
+ k |= cmp( Y, Zero );
+ k |= cmp( Z, Zero );
+ k |= cmp( T, Zero );
+ k |= cmp( StickyBit, Zero );
+ k |= cmp( YY1, Half );
+ if( k == 0 )
+ {
+ RMult = Rounded;
+ printf("Multiplication appears to round correctly.\n");
+ }
+ else
+ {
+ add( X, U2, t );
+ k = cmp( t, Zero );
+ if( cmp( Y, Zero ) >= 0 )
+ k |= 1;
+ add( Z, U2, t );
+ k |= cmp( t, Zero );
+ if( cmp( T, Zero ) >= 0 )
+ k |= 1;
+ add( StickyBit, U2, t );
+ k |= cmp( t, Zero );
+ if( cmp(YY1, Half) >= 0 )
+ k |= 1;
+ if( k == 0 )
+ {
+ printf("Multiplication appears to chop.\n");
+ }
+ else
+ {
+ printf("* is neither chopped nor correctly rounded.\n");
+ }
+ if( (RMult == Rounded) && (GMult == No) )
+ printf("Multiplication has inconsistent result");
+ }
+ }
+else
+ printf("* is neither chopped nor correctly rounded.\n");
+
+/*=============================================*/
+Milestone = 45;
+/*=============================================*/
+add( One, U2, Y2 );
+sub( U2, One, YY1 );
+add( OneAndHalf, U2, Z );
+add( Z, U2, Z );
+div( Y2, Z, X );
+sub( U2, OneAndHalf, T );
+sub( U2, T, T );
+sub( U2, T, Y );
+div( YY1, Y, Y );
+add( Z, U2, Z );
+div( Y2, Z, Z );
+sub( OneAndHalf, X, X );
+sub( T, Y, Y );
+div( YY1, T, T );
+add( OneAndHalf, U2, t );
+sub( t, Z, Z );
+sub( OneAndHalf, U2, t );
+add( t, T, T );
+k = 0;
+if( cmp( X, Zero ) > 0 )
+ k = 1;
+if( cmp( Y, Zero ) > 0 )
+ k = 1;
+if( cmp( Z, Zero ) > 0 )
+ k = 1;
+if( cmp( T, Zero ) > 0 )
+ k = 1;
+if( k == 0 )
+ {
+ div( Y2, OneAndHalf, X );
+ sub( U2, OneAndHalf, Y );
+ add( U2, OneAndHalf, Z );
+ sub( Y, X, X );
+ div( YY1, OneAndHalf, T );
+ div( YY1, Y, Y );
+ add( Z, U2, t );
+ sub( t, T, T );
+ sub( Z, Y, Y );
+ div( Y2, Z, Z );
+ add( Y2, U2, YY1 );
+ div( Y2, YY1, YY1 );
+ sub( OneAndHalf, Z, Z );
+ sub( Y2, YY1, Y2 );
+ sub( U1, F9, YY1 );
+ div( F9, YY1, YY1 );
+ k = cmp( X, Zero );
+ k |= cmp( Y, Zero );
+ k |= cmp( Z, Zero );
+ k |= cmp( T, Zero );
+ k |= cmp( Y2, Zero );
+ sub( Half, YY1, t );
+ sub( Half, F9, t2 );
+ k |= cmp( t, t2 );
+ if( k == 0 )
+ {
+ RDiv = Rounded;
+ printf("Division appears to round correctly.\n");
+ if(GDiv == No)
+ printf("Division test inconsistent\n");
+ }
+ else
+ {
+ k = 0;
+ if( cmp( X, Zero ) >= 0 )
+ k = 1;
+ if( cmp( Y, Zero ) >= 0 )
+ k = 1;
+ if( cmp( Z, Zero ) >= 0 )
+ k = 1;
+ if( cmp( T, Zero ) >= 0 )
+ k = 1;
+ if( cmp( Y, Zero ) >= 0 )
+ k = 1;
+ sub( Half, YY1, t );
+ sub( Half, F9, t2 );
+ if( cmp( t, t2 ) >= 0 )
+ k = 1;
+ if( k == 0 )
+ {
+ RDiv = Chopped;
+ printf("Division appears to chop.\n");
+ }
+ }
+ }
+if(RDiv == Other)
+ printf("/ is neither chopped nor correctly rounded.\n");
+div( Radix, One, BInvrse );
+mul( BInvrse, Radix, t );
+sub( Half, t, t );
+if( cmp( t, Half ) != 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "Radix * ( 1 / Radix ) differs from 1\n" );
+ }
+
+Milestone = 50;
+/*=============================================*/
+add( F9, U1, t );
+sub( Half, t, t );
+k = cmp( t, Half );
+add( BMinusU2, U2, t );
+sub( One, t, t );
+sub( One, Radix, t2 );
+k |= cmp( t, t2 );
+if( k != 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "Incomplete carry-propagation in Addition\n" );
+ }
+mul( U1, U1, X );
+sub( X, One, X );
+sub( U2, One, Y );
+mul( U2, Y, Y );
+add( One, Y, Y );
+sub( Half, F9, Z );
+sub( Half, X, X );
+sub( Z, X, X );
+sub( One, Y, Y );
+if( (cmp(X,Zero) == 0) && (cmp(Y,Zero) == 0) )
+ {
+ RAddSub = Chopped;
+ printf("Add/Subtract appears to be chopped.\n");
+ }
+if(GAddSub == Yes)
+ {
+ add( Half, U2, X );
+ mul( X, U2, X );
+ sub( U2, Half, Y );
+ mul( Y, U2, Y );
+ add( One, X, X );
+ add( One, Y, Y );
+ add( One, U2, t );
+ sub( X, t, X );
+ sub( Y, One, Y );
+ k = cmp(X,Zero);
+ if( k )
+ printf( "1+U2-[u2(1/2+U2)+1] != 0\n" );
+ k2 = cmp(Y,Zero);
+ if( k2 )
+ printf( "1-[U2(1/2-U2)+1] != 0\n" );
+ k |= k2;
+ if( k == 0 )
+ {
+ add( Half, U2, X );
+ mul( X, U1, X );
+ sub( U2, Half, Y );
+ mul( Y, U1, Y );
+ sub( X, One, X );
+ sub( Y, One, Y );
+ sub( X, F9, X );
+ sub( Y, One, Y );
+ k = cmp(X,Zero);
+ if( k )
+ printf( "F9-[1-U1(1/2+U2)] != 0\n" );
+ k2 = cmp(Y,Zero);
+ if( k2 )
+ printf( "1-[1-U1(1/2-U2)] != 0\n" );
+ k |= k2;
+ if( k == 0 )
+ {
+ RAddSub = Rounded;
+ printf("Addition/Subtraction appears to round correctly.\n");
+ if(GAddSub == No)
+ printf( "Add/Subtract test inconsistent\n");
+ }
+ else
+ {
+ printf("Addition/Subtraction neither rounds nor chops.\n");
+ }
+ }
+ else
+ printf("Addition/Subtraction neither rounds nor chops.\n");
+ }
+else
+ printf("Addition/Subtraction neither rounds nor chops.\n");
+
+mov( One, S );
+add( One, Half, X );
+mul( Half, X, X );
+add( One, X, X );
+add( One, U2, Y );
+mul( Y, Half, Y );
+sub( Y, X, Z );
+sub( X, Y, T );
+add( Z, T, StickyBit );
+if( cmp(StickyBit, Zero) != 0 )
+ {
+ mov( Zero, S );
+ ErrCnt[Flaw] += 1;
+ printf( "(X - Y) + (Y - X) is non zero!\n" );
+ }
+mov( Zero, StickyBit );
+FLOOR( RadixD2, t );
+k2 = cmp( t, RadixD2 );
+k = 1;
+if( (GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
+ && (RMult == Rounded) && (RDiv == Rounded)
+ && (RAddSub == Rounded) && (k2 == 0) )
+ {
+ printf("Checking for sticky bit.\n");
+ k = 0;
+ add( Half, U1, X );
+ mul( X, U2, X );
+ mul( Half, U2, Y );
+ add( One, Y, Z );
+ add( One, X, T );
+ sub( One, Z, t );
+ sub( One, T, t2 );
+ if( cmp(t,Zero) > 0 )
+ {
+ k = 1;
+ printf( "[1+(1/2)U2]-1 > 0\n" );
+ }
+ if( cmp(t2,U2) < 0 )
+ {
+ k = 1;
+ printf( "[1+U2(1/2+U1)]-1 < U2\n" );
+ }
+ add( T, Y, Z );
+ sub( X, Z, Y );
+ sub( T, Z, t );
+ sub( T, Y, t2 );
+ if( cmp(t,U2) < 0 )
+ {
+ k = 1;
+ printf( "[[1+U2(1/2+U1)]+(1/2)U2]-[1+U2(1/2+U1)] < U2\n" );
+ }
+ if( cmp(t2,Zero) != 0 )
+ {
+ k = 1;
+ printf( "(1/2)U2-[1+U2(1/2+U1)] != 0\n" );
+ }
+ add( Half, U1, X );
+ mul( X, U1, X );
+ mul( Half, U1, Y );
+ sub( Y, One, Z );
+ sub( X, One, T );
+ sub( One, Z, t );
+ sub( F9, T, t2 );
+ if( cmp(t,Zero) != 0 )
+ {
+ k = 1;
+ printf( "(1-(1/2)U1)-1 != 0\n" );
+ }
+ if( cmp(t2,Zero) != 0 )
+ {
+ k = 1;
+ printf( "[1-U1(1/2+U1)]-F9 != 0\n" );
+ }
+ sub( U1, Half, Z );
+ mul( Z, U1, Z );
+ sub( Z, F9, T );
+ sub( Y, F9, Q );
+ sub( F9, T, t );
+ if( cmp( t, Zero ) != 0 )
+ {
+ k = 1;
+ printf( "[F9-U1(1/2-U1)]-F9 != 0\n" );
+ }
+ sub( U1, F9, t );
+ sub( Q, t, t );
+ if( cmp( t, Zero ) != 0 )
+ {
+ k = 1;
+ printf( "(F9-U1)-(F9-(1/2)U1) != 0\n" );
+ }
+ add( One, U2, Z );
+ mul( Z, OneAndHalf, Z );
+ add( OneAndHalf, U2, T );
+ sub( Z, T, T );
+ add( U2, T, T );
+ div( Radix, Half, X );
+ add( One, X, X );
+ mul( Radix, U2, Y );
+ add( One, Y, Y );
+ mul( X, Y, Z );
+ if( cmp( T, Zero ) != 0 )
+ {
+ k = 1;
+ printf( "(3/2+U2)-3/2(1+U2)+U2 != 0\n" );
+ }
+ mul( Radix, U2, t );
+ add( X, t, t );
+ sub( Z, t, t );
+ if( cmp( t, Zero ) != 0 )
+ {
+ k = 1;
+ printf( "(1+1/2Radix)+Radix*U2-[1+1/(2Radix)][1+Radix*U2] != 0\n" );
+ }
+ if( cmp(Radix, Two) != 0 )
+ {
+ add( Two, U2, X );
+ div( Two, X, Y );
+ sub( One, Y, t );
+ if( cmp( t, Zero) != 0 )
+ k = 1;
+ }
+ }
+if( k == 0 )
+ {
+ printf("Sticky bit apparently used correctly.\n");
+ mov( One, StickyBit );
+ }
+else
+ {
+ printf("Sticky bit used incorrectly or not at all.\n");
+ }
+
+if( GMult == No || GDiv == No || GAddSub == No ||
+ RMult == Other || RDiv == Other || RAddSub == Other)
+ {
+ ErrCnt[Flaw] += 1;
+ printf("lack(s) of guard digits or failure(s) to correctly round or chop\n");
+printf( "(noted above) count as one flaw in the final tally below\n" );
+ }
+/*=============================================*/
+Milestone = 60;
+/*=============================================*/
+printf("\n");
+printf("Does Multiplication commute? ");
+printf("Testing on %d random pairs.\n", NoTrials);
+SQRT( Three, Random9 );
+mov( Third, Random1 );
+I = 1;
+do
+ {
+ Random();
+ mov( Random1, X );
+ Random();
+ mov( Random1, Y );
+ mul( Y, X, Z9 );
+ mul( X, Y, Z );
+ sub( Z9, Z, Z9 );
+ I = I + 1;
+ }
+while ( ! ((I > NoTrials) || (cmp(Z9,Zero) != 0)));
+if(I == NoTrials)
+ {
+ div( Three, Half, t );
+ add( One, t, Random1 );
+ add( U2, U1, t );
+ add( t, One, Random2 );
+ mul( Random1, Random2, Z );
+ mul( Random2, Random1, Y );
+/* Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
+ * Three) * ((U2 + U1) + One);
+ */
+ div( Three, Half, t2 );
+ add( One, t2, t2 );
+ add( U2, U1, t );
+ add( t, One, t );
+ mul( t2, t, Z9 );
+ mul( t2, t, t );
+ sub( t, Z9, Z9 );
+ }
+if(! ((I == NoTrials) || (cmp(Z9,Zero) == 0)))
+ {
+ ErrCnt[Defect] += 1;
+ printf( "X * Y == Y * X trial fails.\n");
+ }
+else
+ {
+ printf(" No failures found in %d integer pairs.\n", NoTrials);
+ }
+/*=============================================*/
+Milestone = 70;
+/*=============================================*/
+sqtest();
+Milestone = 90;
+pow1test();
+
+Milestone = 110;
+
+printf("Seeking Underflow thresholds UfThold and E0.\n");
+mov( U1, D );
+FLOOR( Precision, t );
+if( cmp(Precision, t) != 0 )
+ {
+ mov( BInvrse, D );
+ mov( Precision, X );
+ do
+ {
+ mul( D, BInvrse, D );
+ sub( One, X, X );
+ }
+ while( cmp(X, Zero) > 0 );
+ }
+mov( One, Y );
+mov( D, Z );
+/* ... D is power of 1/Radix < 1. */
+sigsave = sigfpe;
+if( setjmp(ovfl_buf) )
+ goto under0;
+do
+ {
+ mov( Y, C );
+ mov( Z, Y );
+ mul( Y, Y, Z );
+ add( Z, Z, t );
+ }
+while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) );
+
+under0:
+sigsave = 0;
+
+mov( C, Y );
+mul( Y, D, Z );
+sigsave = sigfpe;
+if( setjmp(ovfl_buf) )
+ goto under1;
+do
+ {
+ mov( Y, C );
+ mov( Z, Y );
+ mul( Y, D, Z );
+ add( Z, Z, t );
+ }
+while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) );
+
+under1:
+sigsave = 0;
+
+if( cmp(Radix,Two) < 0 )
+ mov( Two, HInvrse );
+else
+ mov( Radix, HInvrse );
+div( HInvrse, One, H );
+/* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
+div( C, One, CInvrse );
+mov( C, E0 );
+mul( E0, H, Z );
+/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
+sigsave = sigfpe;
+if( setjmp(ovfl_buf) )
+ goto under2;
+do
+ {
+ mov( E0, Y );
+ mov( Z, E0 );
+ mul( E0, H, Z );
+ add( Z, Z, t );
+ }
+while( (cmp(E0,Z) > 0) && (cmp(t,Z) > 0) );
+
+under2:
+sigsave = 0;
+
+mov( E0, UfThold );
+mov( Zero, E1 );
+mov( Zero, Q );
+mov( U2, E9 );
+add( One, E9, S );
+mul( C, S, D );
+if( cmp(D,C) <= 0 )
+ {
+ mul( Radix, U2, E9 );
+ add( One, E9, S );
+ mul( C, S, D );
+ if( cmp(D, C) <= 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "multiplication gets too many last digits wrong.\n" );
+ mov( E0, Underflow );
+ mov( Zero, YY1 );
+ mov( Z, PseudoZero );
+ }
+ }
+else
+ {
+ mov( D, Underflow );
+ mul( Underflow, H, PseudoZero );
+ mov( Zero, UfThold );
+ do
+ {
+ mov( Underflow, YY1 );
+ mov( PseudoZero, Underflow );
+ add( E1, E1, t );
+ if( cmp(t, E1) <= 0)
+ {
+ mul( Underflow, HInvrse, Y2 );
+ sub( Y2, YY1, E1 );
+ FABS( E1 );
+ mov( YY1, Q );
+ if( (cmp( UfThold, Zero ) == 0)
+ && (cmp(YY1, Y2) != 0) )
+ mov( YY1, UfThold );
+ }
+ mul( PseudoZero, H, PseudoZero );
+ add( PseudoZero, PseudoZero, t );
+ }
+ while( (cmp(Underflow, PseudoZero) > 0)
+ && (cmp(t, PseudoZero) > 0) );
+ }
+/* Comment line 4530 .. 4560 */
+if( cmp(PseudoZero, Zero) != 0 )
+ {
+ printf("\n");
+ mov(PseudoZero, Z );
+/* ... Test PseudoZero for "phoney- zero" violates */
+/* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
+ ... */
+ if( cmp(PseudoZero, Zero) <= 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf("Positive expressions can underflow to an\n");
+ printf("allegedly negative value\n");
+ printf("PseudoZero that prints out as: " );
+ show( PseudoZero );
+ mov( PseudoZero, X );
+ neg( X );
+ if( cmp(X, Zero) <= 0 )
+ {
+ printf("But -PseudoZero, which should be\n");
+ printf("positive, isn't; it prints out as " );
+ show( X );
+ }
+ }
+ else
+ {
+ ErrCnt[Flaw] += 1;
+ printf( "Underflow can stick at an allegedly positive\n");
+ printf("value PseudoZero that prints out as " );
+ show( PseudoZero );
+ }
+/* TstPtUf();*/
+ }
+
+/*=============================================*/
+Milestone = 120;
+/*=============================================*/
+mul( CInvrse, Y, t );
+mul( CInvrse, YY1, t2 );
+if( cmp(t,t2) > 0 )
+ {
+ mul( H, S, S );
+ mov( Underflow, E0 );
+ }
+if(! ((cmp(E1,Zero) == 0) || (cmp(E1,E0) == 0)) )
+ {
+ ErrCnt[Defect] += 1;
+ if( cmp(E1,E0) < 0 )
+ {
+ printf("Products underflow at a higher");
+ printf(" threshold than differences.\n");
+ if( cmp(PseudoZero,Zero) == 0 )
+ mov( E1, E0 );
+ }
+ else
+ {
+ printf("Difference underflows at a higher");
+ printf(" threshold than products.\n");
+ }
+ }
+printf("Smallest strictly positive number found is E0 = " );
+show( E0 );
+mov( E0, Z );
+TstPtUf();
+mov( E0, Underflow );
+if(N == 1)
+ mov( Y, Underflow );
+I = 4;
+if( cmp(E1,Zero) == 0 )
+ I = 3;
+if( cmp( UfThold,Zero) == 0 )
+ I = I - 2;
+UfNGrad = True;
+switch(I)
+ {
+ case 1:
+ mov( Underflow, UfThold );
+ mul( CInvrse, Q, t );
+ mul( CInvrse, Y, t2 );
+ mul( t2, S, t2 );
+ if( cmp( t, t2 ) != 0 )
+ {
+ mov( Y, UfThold );
+ ErrCnt[Failure] += 1;
+ printf( "Either accuracy deteriorates as numbers\n");
+ printf("approach a threshold = " );
+ show( UfThold );
+ printf(" coming down from " );
+ show( C );
+ printf(" or else multiplication gets too many last digits wrong.\n");
+ }
+ break;
+
+ case 2:
+ ErrCnt[Failure] += 1;
+ printf( "Underflow confuses Comparison which alleges that\n");
+ printf("Q == Y while denying that |Q - Y| == 0; these values\n");
+ printf("print out as Q = " );
+ show( Q );
+ printf( ", Y = " );
+ show( Y );
+ sub( Y2, Q, t );
+ FABS(t);
+ printf ("|Q - Y| = " );
+ show( t );
+ mov( Q, UfThold );
+ break;
+
+ case 3:
+ mov( X, X );
+ break;
+
+ case 4:
+ div( E9, E1, t );
+ sub( t, UfThold, t );
+ FABS(t);
+ if( (cmp(Q,UfThold) == 0) && (cmp(E1,E0) == 0)
+ && (cmp(t,E1) <= 0) )
+ {
+ UfNGrad = False;
+ printf("Underflow is gradual; it incurs Absolute Error =\n");
+ printf("(roundoff in UfThold) < E0.\n");
+ mul( E0, CInvrse, Y );
+ add( OneAndHalf, U2, t );
+ mul( Y, t, Y );
+ add( One, U2, X );
+ mul( CInvrse, X, X );
+ div( X, Y, t );
+ IEEE = (cmp(t,E0) == 0);
+ if( IEEE == 0 )
+ {
+ printf( "((CInvrse E0) (1.5+U2)) / (CInvrse (1+U2)) != E0\n" );
+ printf( "CInvrse = " );
+ show( CInvrse );
+ printf( "E0 = " );
+ show( E0 );
+ printf( "U2 = " );
+ show( U2 );
+ printf( "X = " );
+ show(X);
+ printf( "Y = " );
+ show(Y);
+ printf( "Y/X = " );
+ show(t);
+ }
+ }
+ }
+if(UfNGrad)
+ {
+ printf("\n");
+ div( UfThold, Underflow, R );
+ SQRT( R, R );
+ if( cmp(R,H) <= 0)
+ {
+ mul( R, UfThold, Z );
+/* X = Z * (One + R * H * (One + H));*/
+ add( One, H, X );
+ mul( H, X, X );
+ mul( R, X, X );
+ add( One, X, X );
+ mul( Z, X, X );
+ }
+ else
+ {
+ mov( UfThold, Z );
+/*X = Z * (One + H * H * (One + H));*/
+ add( One, H, X );
+ mul( H, X, X );
+ mul( H, X, X );
+ add( One, X, X );
+ mul( Z, X, X );
+ }
+ sub( Z, X, t );
+/* if(! ((cmp(X,Z) == 0) || (cmp(t,Zero) != 0)) )*/
+ if( (cmp(X,Z) != 0) && (cmp(t,Zero) == 0) )
+ {
+/* ErrCnt[Flaw] += 1;*/
+ ErrCnt[Serious] += 1;
+ printf("X = " );
+ show( X );
+ printf( "\tis not equal to Z = " );
+ show( Z );
+/* sub( Z, X, Z9 );*/
+ printf("yet X - Z yields " );
+ show( t );
+ printf("which compares equal to " );
+ show( Zero );
+ printf(" Should this NOT signal Underflow, ");
+ printf("this is a SERIOUS DEFECT\nthat causes ");
+ printf("confusion when innocent statements like\n");;
+ printf(" if (X == Z) ... else");
+ printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
+ printf("encounter Division by Zero although actually\n");
+ printf("X / Z = 1 + " );
+ div( Z, X, t );
+ sub( Half, t, t );
+ sub( Half, t, t );
+ show(t);
+ }
+ }
+printf("The Underflow threshold is " );
+show( UfThold );
+printf( "below which calculation may suffer larger Relative error than" );
+printf( " merely roundoff.\n");
+mul( U1, U1, Y2 );
+mul( Y2, Y2, Y );
+mul( Y, U1, Y2 );
+if( cmp( Y2,UfThold) <= 0 )
+ {
+ if( cmp(Y,E0) > 0 )
+ {
+ ErrCnt[Defect] += 1;
+ I = 5;
+ }
+ else
+ {
+ ErrCnt[Serious] += 1;
+ I = 4;
+ }
+ printf("Range is too narrow; U1^%d Underflows.\n", I);
+ }
+Milestone = 130;
+
+/*Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;*/
+LOG( UfThold, Y );
+LOG( HInvrse, t );
+div( t, Y, Y );
+mul( TwoForty, Y, Y );
+sub( Y, Half, Y );
+FLOOR( Y, Y );
+div( TwoForty, Y, Y );
+neg(Y);
+sub( One, Y, Y2 ); /* ***** changed from Y2 = Y + Y */
+printf("Since underflow occurs below the threshold\n");
+printf("UfThold = " );
+show( HInvrse );
+printf( "\tto the power " );
+show( Y );
+printf( "only underflow should afflict the expression " );
+show( HInvrse );
+printf( "\tto the power " );
+show( Y2 );
+POW( HInvrse, Y2, V9 );
+printf("Actually calculating yields: " );
+show( V9 );
+add( Radix, Radix, t );
+add( t, E9, t );
+mul( t, UfThold, t );
+if( (cmp(V9,Zero) < 0) || (cmp(V9,t) > 0) )
+ {
+ ErrCnt[Serious] += 1;
+ printf( "this is not between 0 and underflow\n");
+ printf(" threshold = " );
+ show( UfThold );
+ }
+else
+ {
+ add( One, E9, t );
+ mul( UfThold, t, t );
+ if( cmp(V9,t) <= 0 )
+ printf("This computed value is O.K.\n");
+ else
+ {
+ ErrCnt[Defect] += 1;
+ printf( "this is not between 0 and underflow\n");
+ printf(" threshold = " );
+ show( UfThold );
+ }
+ }
+
+Milestone = 140;
+
+pow2test();
+
+/*=============================================*/
+Milestone = 160;
+/*=============================================*/
+Pause();
+printf("Searching for Overflow threshold:\n");
+printf("This may generate an error.\n");
+sigsave = sigfpe;
+I = 0;
+mov( CInvrse, Y ); /* a large power of 2 */
+neg(Y);
+mul( HInvrse, Y, V9 ); /* HInvrse = 2 */
+if (setjmp(ovfl_buf))
+ goto overflow;
+do
+ {
+ mov( Y, V );
+ mov( V9, Y );
+ mul( HInvrse, Y, V9 );
+ }
+while( cmp(V9,Y) < 0 ); /* V9 = 2 * Y */
+I = 1;
+
+overflow:
+
+show( HInvrse );
+printf( "\ttimes " );
+show( Y );
+printf( "\tequals " );
+show( V9 );
+
+mov( V9, Z );
+printf("Can `Z = -Y' overflow?\n");
+printf("Trying it on Y = " );
+show(Y);
+mov( Y, V9 );
+neg( V9 );
+mov( V9, V0 );
+sub( Y, V, t );
+add( V, V0, t2 );
+if( cmp(t,t2) == 0 )
+ printf("Seems O.K.\n");
+else
+ {
+ printf("finds a Flaw, -(-Y) differs from Y.\n");
+ printf( "V-Y=t:" );
+ show(V);
+ show(Y);
+ show(t);
+ printf( "V+V0=t2:" );
+ show(V);
+ show(V0);
+ show(t2);
+ ErrCnt[Flaw] += 1;
+ }
+if( (cmp(Z, Y) != 0) && (I != 0) )
+ {
+ ErrCnt[Serious] += 1;
+ printf("overflow past " );
+ show( Y );
+ printf( "\tshrinks to " );
+ show( Z );
+ printf( "= Y * " );
+ show( HInvrse );
+ }
+/*Y = V * (HInvrse * U2 - HInvrse);*/
+mul( HInvrse, U2, Y );
+sub( HInvrse, Y, Y );
+mul( V, Y, Y );
+/*Z = Y + ((One - HInvrse) * U2) * V;*/
+sub( HInvrse, One, Z );
+mul( Z, U2, Z );
+mul( Z, V, Z );
+add( Y, Z, Z );
+if( cmp(Z,V0) < 0 )
+ mov( Z, Y );
+if( cmp(Y,V0) < 0)
+ mov( Y, V );
+sub( V, V0, t );
+if( cmp(t,V0) < 0 )
+ mov( V0, V );
+printf("Overflow threshold is V = " );
+show( V );
+if(I)
+ {
+ printf("Overflow saturates at V0 = " );
+ show( V0 );
+ }
+else
+printf("There is no saturation value because the system traps on overflow.\n");
+
+mul( V, One, V9 );
+printf("No Overflow should be signaled for V * 1 = " );
+show( V9 );
+div( One, V, V9 );
+ printf(" nor for V / 1 = " );
+ show( V9 );
+ printf("Any overflow signal separating this * from the one\n");
+ printf("above is a DEFECT.\n");
+/*=============================================*/
+Milestone = 170;
+/*=============================================*/
+mov( V, t );
+neg( t );
+k = 0;
+if( cmp(t,V) >= 0 )
+ k = 1;
+mov( V0, t );
+neg( t );
+if( cmp(t,V0) >= 0 )
+ k = 1;
+mov( UfThold, t );
+neg(t);
+if( cmp(t,V) >= 0 )
+ k = 1;
+if( cmp(UfThold,V) >= 0 )
+ k = 1;
+if( k != 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "Comparisons involving +-");
+ show( V );
+ show( V0 );
+ show( UfThold );
+ printf("are confused by Overflow." );
+ }
+/*=============================================*/
+Milestone = 175;
+/*=============================================*/
+printf("\n");
+for(Indx = 1; Indx <= 3; ++Indx) {
+ switch(Indx)
+ {
+ case 1: mov(UfThold, Z); break;
+ case 2: mov( E0, Z); break;
+ case 3: mov(PseudoZero, Z); break;
+ }
+if( cmp(Z, Zero) != 0 )
+ {
+ SQRT( Z, V9 );
+ mul( V9, V9, Y );
+ mul( Radix, E9, t );
+ sub( t, One, t );
+ div( t, Y, t );
+ add( One, Radix, t2 );
+ add( t2, E9, t2 );
+ mul( t2, Z, t2 );
+ if( (cmp(t,Z) < 0) || (cmp(Y,t2) > 0) )
+ {
+ if( cmp(V9,U1) > 0 )
+ ErrCnt[Serious] += 1;
+ else
+ ErrCnt[Defect] += 1;
+ printf("Comparison alleges that what prints as Z = " );
+ show( Z );
+ printf(" is too far from sqrt(Z) ^ 2 = " );
+ show( Y );
+ }
+ }
+}
+
+Milestone = 180;
+
+for(Indx = 1; Indx <= 2; ++Indx)
+ {
+ if(Indx == 1)
+ mov( V, Z );
+ else
+ mov( V0, Z );
+ SQRT( Z, V9 );
+ mul( Radix, E9, X );
+ sub( X, One, X );
+ mul( X, V9, X );
+ mul( V9, X, V9 );
+ mul( Two, Radix, t );
+ mul( t, E9, t );
+ sub( t, One, t );
+ mul( t, Z, t );
+ if( (cmp(V9,t) < 0) || (cmp(V9,Z) > 0) )
+ {
+ mov( V9, Y );
+ if( cmp(X,W) < 0 )
+ ErrCnt[Serious] += 1;
+ else
+ ErrCnt[Defect] += 1;
+ printf("Comparison alleges that Z = " );
+ show( Z );
+ printf(" is too far from sqrt(Z) ^ 2 :" );
+ show( Y );
+ }
+ }
+
+Milestone = 190;
+
+Pause();
+mul( UfThold, V, X );
+mul( Radix, Radix, Y );
+mul( X, Y, t );
+if( (cmp(t,One) < 0) || (cmp(X,Y) > 0) )
+ {
+ mul( X, Y, t );
+ div( U1, Y, t2 );
+ if( (cmp(t,U1) < 0) || (cmp(X,t2) > 0) )
+ {
+ ErrCnt[Defect] += 1;
+ printf( "Badly " );
+ }
+ else
+ {
+ ErrCnt[Flaw] += 1;
+ }
+ printf(" unbalanced range; UfThold * V = " );
+ show( X );
+ printf( "\tis too far from 1.\n");
+ }
+Milestone = 200;
+
+for(Indx = 1; Indx <= 5; ++Indx)
+ {
+ mov( F9, X );
+ switch(Indx)
+ {
+ case 2: add( One, U2, X ); break;
+ case 3: mov( V, X ); break;
+ case 4: mov(UfThold,X); break;
+ case 5: mov(Radix,X);
+ }
+ mov( X, Y );
+
+ sigsave = sigfpe;
+ if (setjmp(ovfl_buf))
+ {
+ printf(" X / X traps when X = " );
+ show( X );
+ }
+ else
+ {
+/*V9 = (Y / X - Half) - Half;*/
+ div( X, Y, t );
+ sub( Half, t, t );
+ sub( Half, t, V9 );
+ if( cmp(V9,Zero) == 0 )
+ continue;
+ mov( U1, t );
+ neg(t);
+ if( (cmp(V9,t) == 0) && (Indx < 5) )
+ {
+ ErrCnt[Flaw] += 1;
+ }
+ else
+ {
+ ErrCnt[Serious] += 1;
+ }
+ printf(" X / X differs from 1 when X = " );
+ show( X );
+ printf(" instead, X / X - 1/2 - 1/2 = " );
+ show( V9 );
+ }
+ }
+
+ Pause();
+ printf("\n");
+ {
+ static char *msg[] = {
+ "FAILUREs encountered =",
+ "SERIOUS DEFECTs discovered =",
+ "DEFECTs discovered =",
+ "FLAWs discovered =" };
+ int i;
+ for(i = 0; i < 4; i++) if (ErrCnt[i])
+ printf("The number of %-29s %d.\n",
+ msg[i], ErrCnt[i]);
+ }
+ printf("\n");
+ if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
+ + ErrCnt[Flaw]) > 0) {
+ if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
+ Defect] == 0) && (ErrCnt[Flaw] > 0)) {
+ printf("The arithmetic diagnosed seems ");
+ printf("satisfactory though flawed.\n");
+ }
+ if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
+ && ( ErrCnt[Defect] > 0)) {
+ printf("The arithmetic diagnosed may be acceptable\n");
+ printf("despite inconvenient Defects.\n");
+ }
+ if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
+ printf("The arithmetic diagnosed has ");
+ printf("unacceptable serious defects.\n");
+ }
+ if (ErrCnt[Failure] > 0) {
+ printf("Fatal FAILURE may have spoiled this");
+ printf(" program's subsequent diagnoses.\n");
+ }
+ }
+ else {
+ printf("No failures, defects nor flaws have been discovered.\n");
+ if (! ((RMult == Rounded) && (RDiv == Rounded)
+ && (RAddSub == Rounded) && (RSqrt == Rounded)))
+ printf("The arithmetic diagnosed seems satisfactory.\n");
+ else {
+ k = 0;
+ if( cmp( Radix, Two ) == 0 )
+ k = 1;
+ if( cmp( Radix, Ten ) == 0 )
+ k = 1;
+ if( (cmp(StickyBit,One) >= 0) && (k == 1) )
+ {
+ printf("Rounding appears to conform to ");
+ printf("the proposed IEEE standard P");
+ k = 0;
+ k |= cmp( Radix, Two );
+ mul( Four, Three, t );
+ mul( t, Two, t );
+ sub( t, Precision, t );
+ sub( TwentySeven, Precision, t2 );
+ sub( TwentySeven, t2, t2 );
+ add( t2, One, t2 );
+ mul( t2, t, t );
+ if( (cmp(Radix,Two) == 0)
+ && (cmp(t,Zero) == 0) )
+ printf("754");
+ else
+ printf("854");
+ if(IEEE)
+ printf(".\n");
+ else
+ {
+ printf(",\nexcept for possibly Double Rounding");
+ printf(" during Gradual Underflow.\n");
+ }
+ }
+ printf("The arithmetic diagnosed appears to be excellent!\n");
+ }
+ }
+ if (fpecount)
+ printf("\nA total of %d floating point exceptions were registered.\n",
+ fpecount);
+ printf("END OF TEST.\n");
+ }
+
+
+/* Random */
+/* Random computes
+ X = (Random1 + Random9)^5
+ Random1 = X - FLOOR(X) + 0.000005 * X;
+ and returns the new value of Random1
+*/
+
+
+static int randflg = 0;
+FLOAT(C5em6);
+
+Random()
+{
+
+if( randflg == 0 )
+ {
+ mov( Six, t );
+ neg(t);
+ POW( Ten, t, t );
+ mul( Five, t, C5em6 );
+ randflg = 1;
+ }
+add( Random1, Random9, t );
+mul( t, t, t2 );
+mul( t2, t2, t2 );
+mul( t, t2, t );
+FLOOR(t, t2 );
+sub( t2, t, t2 );
+mul( t, C5em6, t );
+add( t, t2, Random1 );
+/*return(Random1);*/
+}
+
+/* SqXMinX */
+
+SqXMinX( ErrKind )
+int ErrKind;
+{
+mul( X, BInvrse, t2 );
+sub( t2, X, t );
+/*SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;*/
+mul( X, X, Sqarg );
+SQRT( Sqarg, SqEr );
+sub( t2, SqEr, SqEr );
+sub( t, SqEr, SqEr );
+div( OneUlp, SqEr, SqEr );
+if( cmp(SqEr,Zero) != 0)
+ {
+ Showsq( 0 );
+ add( J, One, J );
+ ErrCnt[ErrKind] += 1;
+ printf("sqrt of " );
+ mul( X, X, t );
+ show( t );
+ printf( "minus " );
+ show( X );
+ printf( "equals " );
+ mul( OneUlp, SqEr, t );
+ show( t );
+ printf("\tinstead of correct value 0 .\n");
+ }
+}
+
+/* NewD */
+
+NewD()
+{
+mul( Z1, Q, X );
+/*X = FLOOR(Half - X / Radix) * Radix + X;*/
+div( Radix, X, t );
+sub( t, Half, t );
+FLOOR( t, t );
+mul( t, Radix, t );
+add( t, X, X );
+/*Q = (Q - X * Z) / Radix + X * X * (D / Radix);*/
+mul( X, Z, t );
+sub( t, Q, t );
+div( Radix, t, t );
+div( Radix, D, t2 );
+mul( X, t2, t2 );
+mul( X, t2, t2 );
+add( t, t2, Q );
+/*Z = Z - Two * X * D;*/
+mul( Two, X, t );
+mul( t, D, t );
+sub( t, Z, Z );
+
+if( cmp(Z,Zero) <= 0)
+ {
+ neg(Z);
+ neg(Z1);
+ }
+mul( Radix, D, D );
+}
+
+/* SR3750 */
+
+SR3750()
+{
+sub( Radix, X, t );
+sub( Radix, Z2, t2 );
+k = 0;
+if( cmp(t,t2) < 0 )
+ k = 1;
+sub( Z2, X, t );
+sub( Z2, W, t2 );
+if( cmp(t,t2) > 0 )
+ k = 1;
+/*if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {*/
+if( k == 0 )
+ {
+ I = I + 1;
+ mul( X, D, X2 );
+ mov( X2, Sqarg );
+ SQRT( X2, X2 );
+/*Y2 = (X2 - Z2) - (Y - Z2);*/
+ sub( Z2, X2, Y2 );
+ sub( Z2, Y, t );
+ sub( t, Y2, Y2 );
+ sub( Half, Y, X2 );
+ div( X2, X8, X2 );
+ mul( Half, X2, t );
+ mul( t, X2, t );
+ sub( t, X2, X2 );
+/*SqEr = (Y2 + Half) + (Half - X2);*/
+ add( Y2, Half, SqEr );
+ sub( X2, Half, t );
+ add( t, SqEr, SqEr );
+ Showsq( -1 );
+ sub( X2, Y2, SqEr );
+ Showsq( 1 );
+ }
+}
+
+/* IsYeqX */
+
+IsYeqX()
+{
+if( cmp(Y,X) != 0 )
+ {
+ if (N <= 0)
+ {
+ if( (cmp(Z,Zero) == 0) && (cmp(Q,Zero) <= 0) )
+ printf("WARNING: computing\n");
+ else
+ {
+ ErrCnt[Defect] += 1;
+ printf( "computing\n");
+ }
+ show( Z );
+ printf( "\tto the power " );
+ show( Q );
+ printf("\tyielded " );
+ show( Y );
+ printf("\twhich compared unequal to correct " );
+ show( X );
+ sub( X, Y, t );
+ printf("\t\tthey differ by " );
+ show( t );
+ }
+ N = N + 1; /* ... count discrepancies. */
+ }
+}
+
+/* SR3980 */
+
+SR3980()
+{
+long li;
+
+do
+ {
+/*Q = (FLOAT) I;*/
+ li = I;
+ LTOF( &li, Q );
+ POW( Z, Q, Y );
+ IsYeqX();
+ if(++I > M)
+ break;
+ mul( Z, X, X );
+ }
+while( cmp(X,W) < 0 );
+}
+
+/* PrintIfNPositive */
+
+PrintIfNPositive()
+{
+if(N > 0)
+ printf("Similar discrepancies have occurred %d times.\n", N);
+}
+
+
+/* TstPtUf */
+
+TstPtUf()
+{
+N = 0;
+if( cmp(Z,Zero) != 0)
+ {
+ printf( "Z = " );
+ show(Z);
+ printf("Since comparison denies Z = 0, evaluating ");
+ printf("(Z + Z) / Z should be safe.\n");
+ sigsave = sigfpe;
+ if (setjmp(ovfl_buf))
+ goto very_serious;
+ add( Z, Z, Q9 );
+ div( Z, Q9, Q9 );
+ printf("What the machine gets for (Z + Z) / Z is " );
+ show( Q9 );
+ sub( Two, Q9, t );
+ FABS(t);
+ mul( Radix, U2, t2 );
+ if( cmp(t,t2) < 0 )
+ {
+ printf("This is O.K., provided Over/Underflow");
+ printf(" has NOT just been signaled.\n");
+ }
+ else
+ {
+ if( (cmp(Q9,One) < 0) || (cmp(Q9,Two) > 0) )
+ {
+very_serious:
+ N = 1;
+ ErrCnt [Serious] = ErrCnt [Serious] + 1;
+ printf("This is a VERY SERIOUS DEFECT!\n");
+ }
+ else
+ {
+ N = 1;
+ ErrCnt[Defect] += 1;
+ printf("This is a DEFECT!\n");
+ }
+ }
+ mul( Z, One, V9 );
+ mov( V9, Random1 );
+ mul( One, Z, V9 );
+ mov( V9, Random2 );
+ div( One, Z, V9 );
+ if( (cmp(Z,Random1) == 0) && (cmp(Z,Random2) == 0)
+ && (cmp(Z,V9) == 0) )
+ {
+ if (N > 0)
+ Pause();
+ }
+ else
+ {
+ N = 1;
+ ErrCnt[Defect] += 1;
+ printf( "What prints as Z = ");
+ show( Z );
+ printf( "\tcompares different from " );
+ if( cmp(Z,Random1) != 0)
+ {
+ printf("Z * 1 = " );
+ show( Random1 );
+ }
+ if( (cmp(Z,Random2) != 0)
+ || (cmp(Random2,Random1) != 0) )
+ {
+ printf("1 * Z == " );
+ show( Random2 );
+ }
+ if( cmp(Z,V9) != 0 )
+ {
+ printf("Z / 1 = " );
+ show( V9 );
+ }
+ if( cmp(Random2,Random1) != 0 )
+ {
+ ErrCnt[Defect] += 1;
+ printf( "Multiplication does not commute!\n");
+ printf("\tComparison alleges that 1 * Z = " );
+ show(Random2);
+ printf("\tdiffers from Z * 1 = " );
+ show(Random1);
+ }
+ Pause();
+ }
+ }
+}
+
+Pause()
+{
+}
+
+Sign( x, y )
+FSIZE *x, *y;
+{
+
+if( cmp( x, Zero ) < 0 )
+ {
+ mov( One, y );
+ neg( y );
+ }
+else
+ {
+ mov( One, y );
+ }
+}
+
+sqtest()
+{
+printf("\nRunning test of square root(x).\n");
+
+RSqrt = Other;
+k = 0;
+SQRT( Zero, t );
+k |= cmp( Zero, t );
+mov( Zero, t );
+neg(t);
+SQRT( t, t2 );
+k |= cmp( t, t2 );
+SQRT( One, t );
+k |= cmp( One, t );
+if( k != 0 )
+ {
+ ErrCnt[Failure] += 1;
+ printf( "Square root of 0.0, -0.0 or 1.0 wrong\n");
+ }
+mov( Zero, MinSqEr );
+mov( Zero, MaxSqEr );
+mov( Zero, J );
+mov( Radix, X );
+mov( U2, OneUlp );
+SqXMinX( Serious );
+mov( BInvrse, X );
+mul( BInvrse, U1, OneUlp );
+SqXMinX( Serious );
+mov( U1, X );
+mul( U1, U1, OneUlp );
+SqXMinX( Serious );
+if( cmp(J,Zero) != 0)
+ Pause();
+printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
+mov( Zero, J );
+mov( Two, X );
+mov( Radix, Y );
+if( cmp(Radix,One) != 0 )
+ {
+ lngint = NoTrials;
+ LTOF( &lngint, t );
+ FTOL( t, &lng2, X );
+ if( lngint != lng2 )
+ {
+ printf( "Integer conversion error\n" );
+ exit(1);
+ }
+ do
+ {
+ mov( Y, X );
+ mul( Radix, Y, Y );
+ sub( X, Y, t2 );
+ }
+ while( ! (cmp(t2,t) >= 0) );
+ }
+mul( X, U2, OneUlp );
+I = 1;
+while(I < 10)
+ {
+ add( X, One, X );
+ SqXMinX( Defect );
+ if( cmp(J,Zero) > 0 )
+ break;
+ I = I + 1;
+ }
+printf("Test for sqrt monotonicity.\n");
+I = - 1;
+mov( BMinusU2, X );
+mov( Radix, Y );
+mul( Radix, U2, Z );
+add( Radix, Z, Z );
+NotMonot = False;
+Monot = False;
+while( ! (NotMonot || Monot))
+ {
+ I = I + 1;
+ SQRT(X, X);
+ SQRT(Y,Q);
+ SQRT(Z,Z);
+ if( (cmp(X,Q) > 0) || (cmp(Q,Z) > 0) )
+ NotMonot = True;
+ else
+ {
+ add( Q, Half, Q );
+ FLOOR( Q, Q );
+ mul( Q, Q, t );
+ if( (I > 0) || (cmp(Radix,t) == 0) )
+ Monot = True;
+ else if (I > 0)
+ {
+ if(I > 1)
+ Monot = True;
+ else
+ {
+ mul( Y, BInvrse, Y );
+ sub( U1, Y, X );
+ add( Y, U1, Z );
+ }
+ }
+ else
+ {
+ mov( Q, Y );
+ sub( U2, Y, X );
+ add( Y, U2, Z );
+ }
+ }
+ }
+if( Monot )
+ printf("sqrt has passed a test for Monotonicity.\n");
+else
+ {
+ ErrCnt[Defect] += 1;
+ printf("sqrt(X) is non-monotonic for X near " );
+ show(Y);
+ }
+/*=============================================*/
+Milestone = 80;
+/*=============================================*/
+add( MinSqEr, Half, MinSqEr );
+sub( Half, MaxSqEr, MaxSqEr);
+/*Y = (SQRT(One + U2) - One) / U2;*/
+add( One, U2, Sqarg );
+SQRT( Sqarg, Y );
+sub( One, Y, Y );
+div( U2, Y, Y );
+/*SqEr = (Y - One) + U2 / Eight;*/
+sub( One, Y, t );
+div( Eight, U2, SqEr );
+add( t, SqEr, SqEr );
+Showsq( 1 );
+div( Eight, U2, SqEr );
+add( Y, SqEr, SqEr );
+Showsq( -1 );
+/*Y = ((SQRT(F9) - U2) - (One - U2)) / U1;*/
+mov( F9, Sqarg );
+SQRT( Sqarg, Y );
+sub( U2, Y, Y );
+sub( U2, One, t );
+sub( t, Y, Y );
+div( U1, Y, Y );
+div( Eight, U1, SqEr );
+add( Y, SqEr, SqEr );
+Showsq( 1 );
+/*SqEr = (Y + One) + U1 / Eight;*/
+div( Eight, U1, t );
+add( Y, One, SqEr );
+add( SqEr, t, SqEr );
+Showsq( -1 );
+mov( U2, OneUlp );
+mov( OneUlp, X );
+for( Indx = 1; Indx <= 3; ++Indx)
+ {
+/*Y = SQRT((X + U1 + X) + F9);*/
+ add( X, U1, Y );
+ add( Y, X, Y );
+ add( Y, F9, Y );
+ mov( Y, Sqarg );
+ SQRT( Sqarg, Y );
+/*Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;*/
+ sub( U2, One, t );
+ add( t, X, t );
+ sub( U2, Y, Y );
+ sub( t, Y, Y );
+ div( OneUlp, Y, Y );
+/*Z = ((U1 - X) + F9) * Half * X * X / OneUlp;*/
+ sub( X, U1, t );
+ add( t, F9, t );
+ mul( t, Half, t );
+ mul( t, X, t );
+ mul( t, X, t );
+ div( OneUlp, t, Z );
+ add( Y, Half, SqEr );
+ add( SqEr, Z, SqEr );
+ Showsq( -1 );
+ sub( Half, Y, SqEr );
+ add( SqEr, Z, SqEr );
+ Showsq( 1 );
+ if(((Indx == 1) || (Indx == 3)))
+ {
+/*X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));*/
+ mov( OneUlp, Sqarg );
+ SQRT( Sqarg, t );
+ mul( Nine, t, t );
+ div( t, Eight, t );
+ FLOOR( t, t );
+ Sign( X, t2 );
+ mul( t2, t, t );
+ mul( OneUlp, t, X );
+ }
+ else
+ {
+ mov( U1, OneUlp );
+ mov( OneUlp, X );
+ neg( X );
+ }
+ }
+/*=============================================*/
+Milestone = 85;
+/*=============================================*/
+SqRWrng = False;
+Anomaly = False;
+if( cmp(Radix,One) != 0 )
+ {
+ printf("Testing whether sqrt is rounded or chopped.\n");
+/*D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));*/
+ FLOOR( Precision, t2 );
+ add( One, Precision, t );
+ sub( t2, t, t );
+ POW( Radix, t, D );
+ add( Half, D, D );
+ FLOOR( D, D );
+/* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
+ div( Radix, D, X );
+ div( A1, D, Y );
+ FLOOR( X, t );
+ FLOOR( Y, t2 );
+ if( (cmp(X,t) != 0) || (cmp(Y,t2) != 0) )
+ {
+ Anomaly = True;
+ printf( "Anomaly 1\n" );
+ }
+ else
+ {
+ mov( Zero, X );
+ mov( X, Z2 );
+ mov( One, Y );
+ mov( Y, Y2 );
+ sub( One, Radix, Z1 );
+ mul( Four, D, FourD );
+ do
+ {
+ if( cmp(Y2,Z2) >0 )
+ {
+ mov( Radix, Q );
+ mov( Y, YY1 );
+ do
+ {
+/*X1 = FABS(Q + FLOOR(Half - Q / YY1) * YY1);*/
+ div( YY1, Q, t );
+ sub( t, Half, t );
+ FLOOR( t, t );
+ mul( t, YY1, t );
+ add( Q, t, X1 );
+ FABS( X1 );
+ mov( YY1, Q );
+ mov( X1, YY1 );
+ }
+ while( ! (cmp(X1,Zero) <= 0) );
+ if( cmp(Q,One) <= 0 )
+ {
+ mov( Y2, Z2 );
+ mov( Y, Z );
+ }
+ }
+ add( Y, Two, Y );
+ add( X, Eight, X );
+ add( Y2, X, Y2 );
+ if( cmp(Y2,FourD) >= 0 )
+ sub( FourD, Y2, Y2 );
+ }
+ while( ! (cmp(Y,D) >= 0) );
+ sub( Z2, FourD, X8 );
+ mul( Z, Z, Q );
+ add( X8, Q, Q );
+ div( FourD, Q, Q );
+ div( Eight, X8, X8 );
+ FLOOR( Q, t );
+ if( cmp(Q,t) != 0 )
+ {
+ Anomaly = True;
+ printf( "Anomaly 2\n" );
+ }
+ else
+ {
+ Break = False;
+ do
+ {
+ mul( Z1, Z, X );
+/*X = X - FLOOR(X / Radix) * Radix;*/
+ div( Radix, X, t );
+ FLOOR( t, t );
+ mul( t, Radix, t );
+ sub( t, X, X );
+ if( cmp(X,One) == 0 )
+ Break = True;
+ else
+ sub( One, Z1, Z1 );
+ }
+ while( ! (Break || (cmp(Z1,Zero) <= 0)) );
+ if( (cmp(Z1,Zero) <= 0) && (! Break))
+ {
+ printf( "Anomaly 3\n" );
+ Anomaly = True;
+ }
+ else
+ {
+ if( cmp(Z1,RadixD2) > 0)
+ sub( Radix, Z1, Z1 );
+ do
+ {
+ NewD();
+ mul( U2, D, t );
+ }
+ while( ! (cmp(t,F9) >= 0) );
+ mul( D, Radix, t );
+ sub( D, t, t );
+ sub( D, W, t2 );
+ if (cmp(t,t2) != 0 )
+ {
+ printf( "Anomaly 4\n" );
+ Anomaly = True;
+ }
+ else
+ {
+ mov( D, Z2 );
+ I = 0;
+ add( One, Z, t );
+ mul( t, Half, t );
+ add( D, t, Y );
+ add( D, Z, X );
+ add( X, Q, X );
+ SR3750();
+ sub( Z, One, t );
+ mul( t, Half, t );
+ add( D, t, Y );
+ add( Y, D, Y );
+ sub( Z, D, X );
+ add( X, D, X );
+ add( X, Q, t );
+ add( t, X, X );
+ SR3750();
+ NewD();
+ sub( Z2, D, t );
+ sub( Z2, W, t2 );
+ if(cmp(t,t2) != 0 )
+ {
+ printf( "Anomaly 5\n" );
+ Anomaly = True;
+ }
+ else
+ {
+/*Y = (D - Z2) + (Z2 + (One - Z) * Half);*/
+ sub( Z, One, t );
+ mul( t, Half, t );
+ add( Z2, t, t );
+ sub( Z2, D, Y );
+ add( Y, t, Y );
+/*X = (D - Z2) + (Z2 - Z + Q);*/
+ sub( Z, Z2, t );
+ add( t, Q, t );
+ sub( Z2, D, X );
+ add( X, t, X );
+ SR3750();
+ add( One, Z, Y );
+ mul( Y, Half, Y );
+ mov( Q, X );
+ SR3750();
+ if(I == 0)
+ {
+ printf( "Anomaly 6\n" );
+ Anomaly = True;
+ }
+ }
+ }
+ }
+ }
+ }
+ if ((I == 0) || Anomaly)
+ {
+ ErrCnt[Failure] += 1;
+ printf( "Anomalous arithmetic with Integer < \n");
+ printf("Radix^Precision = " );
+ show( W );
+ printf(" fails test whether sqrt rounds or chops.\n");
+ SqRWrng = True;
+ }
+ }
+if(! Anomaly)
+ {
+ if(! ((cmp(MinSqEr,Zero) < 0) || (cmp(MaxSqEr,Zero) > 0))) {
+ RSqrt = Rounded;
+ printf("Square root appears to be correctly rounded.\n");
+ }
+ else
+ {
+ k = 0;
+ add( MaxSqEr, U2, t );
+ sub( Half, U2, t2 );
+ if( cmp(t,t2) > 0 )
+ k = 1;
+ if( cmp( MinSqEr, Half ) > 0 )
+ k = 1;
+ add( MinSqEr, Radix, t );
+ if( cmp( t, Half ) < 0 )
+ k = 1;
+ if( k == 1 )
+ SqRWrng = True;
+ else
+ {
+ RSqrt = Chopped;
+ printf("Square root appears to be chopped.\n");
+ }
+ }
+ }
+if( SqRWrng )
+ {
+ printf("Square root is neither chopped nor correctly rounded.\n");
+ printf("Observed errors run from " );
+ sub( Half, MinSqEr, t );
+ show( t );
+ printf("\tto " );
+ add( Half, MaxSqEr, t );
+ show( t );
+ printf( "ulps.\n" );
+ sub( MinSqEr, MaxSqEr, t );
+ mul( Radix, Radix, t2 );
+ if( cmp( t, t2 ) >= 0 )
+ {
+ ErrCnt[Serious] += 1;
+ printf( "sqrt gets too many last digits wrong\n");
+ }
+ }
+}
+
+Showsq( arg )
+int arg;
+{
+
+k = 0;
+if( arg <= 0 )
+ {
+ if( cmp(SqEr,MinSqEr) < 0 )
+ {
+ k = 1;
+ mov( SqEr, MinSqEr );
+ }
+ }
+if( arg >= 0 )
+ {
+ if( cmp(SqEr,MaxSqEr) > 0 )
+ {
+ k = 2;
+ mov( SqEr, MaxSqEr );
+ }
+ }
+#if DEBUG
+if( k != 0 )
+ {
+ printf( "Square root of " );
+ show( arg );
+ printf( "\tis in error by " );
+ show( SqEr );
+ }
+#endif
+}
+
+
+pow1test()
+{
+
+/*=============================================*/
+Milestone = 90;
+/*=============================================*/
+Pause();
+printf("Testing powers Z^i for small Integers Z and i.\n");
+N = 0;
+/* ... test powers of zero. */
+I = 0;
+mov( Zero, Z );
+neg(Z);
+M = 3;
+Break = False;
+do
+ {
+ mov( One, X );
+ SR3980();
+ if(I <= 10)
+ {
+ I = 1023;
+ SR3980();
+ }
+ if( cmp(Z,MinusOne) == 0 )
+ Break = True;
+ else
+ {
+ mov( MinusOne, Z );
+ PrintIfNPositive();
+ N = 0;
+/* .. if(-1)^N is invalid, replace MinusOne by One. */
+ I = - 4;
+ }
+ }
+while( ! Break );
+PrintIfNPositive();
+N1 = N;
+N = 0;
+mov( A1, Z );
+/*M = FLOOR(Two * LOG(W) / LOG(A1));*/
+LOG( W, t );
+mul( Two, t, t );
+FLOOR( t, t );
+LOG( A1, t2 );
+div( t2, t, t );
+FTOL( t, &lngint, t2 );
+M = lngint;
+Break = False;
+do
+ {
+ mov( Z, X );
+ I = 1;
+ SR3980();
+ if( cmp(Z,AInvrse) == 0 )
+ Break = True;
+ else
+ mov( AInvrse, Z );
+ }
+while( ! (Break) );
+/*=============================================*/
+Milestone = 100;
+/*=============================================*/
+/* Powers of Radix have been tested, */
+/* next try a few primes */
+
+M = NoTrials;
+
+mov( Three, Z );
+do
+ {
+ mov( Z, X );
+ I = 1;
+ SR3980();
+ do
+ {
+ add( Z, Two, Z );
+ div( Three, Z, t );
+ FLOOR( t, t );
+ mul( Three, t, t );
+ }
+ while( cmp(t,Z) == 0 );
+ mul( Eight, Three, t );
+ }
+while( cmp(Z,t) < 0 );
+
+if(N > 0)
+ {
+ printf("Errors like this may invalidate financial calculations\n");
+ printf("\tinvolving interest rates.\n");
+ }
+PrintIfNPositive();
+N += N1;
+if(N == 0)
+ printf("... no discrepancies found.\n");
+if(N > 0)
+ Pause();
+else printf("\n");
+}
+
+
+
+pow2test()
+{
+printf("\n");
+/* ...calculate Exp2 == exp(2) == 7.38905 60989 30650 22723 04275-... */
+mov( Zero, X );
+mov( Two, t2 ); /*I = 2;*/
+
+mul( Two, Three, Y );
+mov( Zero, Q );
+N = 0;
+do
+ {
+ mov( X, Z );
+ add( t2, One, t2 ); /*I = I + 1;*/
+ add( t2, t2, t );
+ div( t, Y, Y ); /*Y = Y / (I + I);*/
+ add( Y, Q, R );
+ add( Z, R, X );
+ sub( X, Z, Q );
+ add( Q, R, Q );
+ }
+while( cmp(X,Z) > 0 );
+
+/*Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);*/
+div( Eight, One, t );
+add( OneAndHalf, t, Z );
+mul( OneAndHalf, ThirtyTwo, t );
+div( t, X, t );
+add( Z, t, Z );
+mul( Z, Z, X );
+mul( X, X, Exp2 );
+mov( F9, X );
+sub( U1, X, Y );
+printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = " );
+show( Exp2 );
+printf( "\tas X -> 1.\n" );
+for(I = 1;;)
+ {
+ sub( BInvrse, X, Z );
+/*Z = (X + One) / (Z - (One - BInvrse));*/
+ add( X, One, t2 );
+ sub( BInvrse, One, t );
+ sub( t, Z, t );
+ div( t, t2, Z );
+ POW( X, Z, Sqarg );
+ sub( Exp2, Sqarg, Q );
+ mov( Q, t );
+ FABS( t );
+ mul( TwoForty, U2, t2 );
+ if( cmp( t, t2 ) > 0 )
+ {
+ N = 1;
+ sub( BInvrse, X, V9 );
+ sub( BInvrse, One, t );
+ sub( t, V9, V9 );
+ ErrCnt[Defect] += 1;
+ printf( "Calculated " );
+ show( Sqarg );
+ printf(" for \t(1 + " );
+ show( V9 );
+ printf( "\tto the power " );
+ show( Z );
+ printf("\tdiffers from correct value by " );
+ show( Q );
+ printf("\tThis much error may spoil financial\n");
+ printf("\tcalculations involving tiny interest rates.\n");
+ break;
+ }
+ else
+ {
+ sub( X, Y, Z );
+ mul( Z, Two, Z );
+ add( Z, Y, Z );
+ mov( Y, X );
+ mov( Z, Y );
+ sub( F9, X, Z );
+ mul( Z, Z, Z );
+ add( Z, One, Z );
+ if( (cmp(Z,One) > 0) && (I < NoTrials) )
+ I++;
+ else
+ {
+ if( cmp(X,One) > 0 )
+ {
+ if(N == 0)
+ printf("Accuracy seems adequate.\n");
+ break;
+ }
+ else
+ {
+ add( One, U2, X );
+ add( U2, U2, Y );
+ add( X, Y, Y );
+ I = 1;
+ }
+ }
+ }
+ }
+/*=============================================*/
+Milestone = 150;
+/*=============================================*/
+printf("Testing powers Z^Q at four nearly extreme values.\n");
+N = 0;
+mov( A1, Z );
+/*Q = FLOOR(Half - LOG(C) / LOG(A1));*/
+LOG( C, t );
+LOG( A1, t2 );
+div( t2, t, t );
+sub( t, Half, t );
+FLOOR( t, Q );
+Break = False;
+do
+ {
+ mov( CInvrse, X );
+ POW( Z, Q, Y );
+ IsYeqX();
+ neg(Q);
+ mov( C, X );
+ POW( Z, Q, Y );
+ IsYeqX();
+ if( cmp(Z,One) < 0 )
+ Break = True;
+ else
+ mov( AInvrse, Z );
+ }
+while( ! (Break));
+PrintIfNPositive();
+if(N == 0)
+ printf(" ... no discrepancies found.\n");
+printf("\n");
+}
|