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