1 /* paranoia.c arithmetic tester
\r
3 * This is an implementation of the PARANOIA program. It substitutes
\r
4 * subroutine calls for ALL floating point arithmetic operations.
\r
5 * This permits you to substitute your own experimental versions of
\r
6 * arithmetic routines. It also defeats compiler optimizations,
\r
7 * so for native arithmetic you can be pretty sure you are testing
\r
8 * the arithmetic and not the compiler.
\r
10 * This version of PARANOIA omits the display of division by zero.
\r
11 * It also omits the test for extra precise subexpressions, since
\r
12 * they cannot occur in this context. Otherwise it includes all the
\r
13 * tests of the 27 Jan 86 distribution, plus a few additional tests.
\r
14 * Commentary has been reduced to a minimum in order to make the program
\r
17 * The original PARANOIA program, written by W. Kahan, C version
\r
18 * by Thos Sumner and David Gay, can be downloaded free from the
\r
19 * Internet NETLIB. An MSDOS disk can be obtained for $15 from:
\r
21 * 6521 Raymond Street
\r
24 * Steve Moshier, 28 Oct 88
\r
25 * last rev: 23 May 92
\r
30 /* To use the native arithmetic of the computer, define NATIVE
\r
31 * to be 1. To use your own supplied arithmetic routines, NATIVE is 0.
\r
35 /* gcc real.c interface */
\r
36 #define L128DOUBLE 0
\r
43 /* Data structure of floating point number. If NATIVE was
\r
44 * selected above, you can define LDOUBLE 1 to test 80-bit long double
\r
45 * precision or define it 0 to test 64-bit double precision.
\r
52 #define FSIZE long double
\r
53 #define FLOAT(x) FSIZE x[NE]
\r
54 static FSIZE eone[NE] = {1.0L}; /* The constant 1.0 */
\r
57 #define ZFLOOR floorl
\r
59 long double sqrtl(), logl(), floorl(), powl();
\r
60 #define FSETUP einit
\r
61 #else /* not LDOUBLE */
\r
62 #define FSIZE double
\r
63 #define FLOAT(x) FSIZE x[NE]
\r
64 static FSIZE eone[NE] = {1.0}; /* The constant 1.0 */
\r
67 #define ZFLOOR floor
\r
69 double sqrt(), log(), floor(), pow();
\r
70 /* Coprocessor initialization,
\r
71 * defeat underflow trap or what have you.
\r
72 * This is required mainly on i386 and 68K processors.
\r
74 #define FSETUP dprec
\r
75 #endif /* double, not LDOUBLE */
\r
77 #else /* not NATIVE */
\r
79 /* Setup for extended double type.
\r
80 * Put NE = 10 for real.c operating with TFmode support (16-byte reals)
\r
81 * Put NE = 6 for real.c operating with XFmode support (10- or 12-byte reals)
\r
82 * The value of NE must agree with that in ehead.h, if ieee.c is used.
\r
85 #define FSIZE unsigned short
\r
86 #define FLOAT(x) unsigned short x[NE]
\r
87 extern unsigned short eone[];
\r
88 #define FSETUP einit
\r
90 /* default for FSETUP */
\r
99 printf( "error: %s\n", s );
\r
102 #endif /* not NATIVE */
\r
107 /* real.c interface */
\r
110 #define FSETUP efsetup
\r
116 /* Use emov to convert from widest type to widest type, ... */
\r
122 /* ... else choose e24toe, e53toe, etc. */
\r
123 #define ENTOE e64toe
\r
124 #define ETOEN etoe64
\r
127 #define NIBITS ((NE-1)*16)
\r
133 ETOEN(eone, enone);
\r
141 unsigned short aa[10], bb[10], cc[10];
\r
154 unsigned short aa[10], bb[10], cc[10];
\r
167 unsigned short aa[10], bb[10], cc[10];
\r
180 unsigned short aa[10], bb[10], cc[10];
\r
192 unsigned short aa[10], bb[10];
\r
208 for( i=0; i<NE; i++ )
\r
216 unsigned short aa[10];
\r
228 for( i=0; i<NE; i++ )
\r
235 unsigned short aa[10];
\r
246 unsigned short aa[10], bb[10];
\r
257 unsigned short aa[10], bb[10];
\r
273 unsigned short aa[10], bb[10], cc[10];
\r
289 unsigned short aa[10], bb[10];
\r
301 unsigned short xx[10], ff[10];
\r
312 unsigned short xx[10];
\r
322 unsigned short xx[10];
\r
328 #else /* not L128DOUBLE */
\r
332 /* Note all arguments of operation subroutines are pointers. */
\r
334 #define add(a,b,c) eadd(a,b,c)
\r
336 #define sub(a,b,c) esub(a,b,c)
\r
338 #define mul(a,b,c) emul(a,b,c)
\r
340 #define div(a,b,c) ediv(a,b,c)
\r
341 /* 1 if a>b, 0 if a==b, -1 if a<b */
\r
342 #define cmp(a,b) ecmp(a,b)
\r
344 #define mov(a,b) emov(a,b)
\r
346 #define neg(a) eneg(a)
\r
348 #define clear(a) eclear(a)
\r
350 #define FABS(x) eabs(x)
\r
351 #define FLOOR(x,y) efloor(x,y)
\r
352 #define LOG(x,y) elog(x,y)
\r
353 #define POW(x,y,z) epow(x,y,z)
\r
354 #define SQRT(x,y) esqrt(x,y)
\r
356 /* x = &FLOAT input, i = &long integer part, f = &FLOAT fractional part */
\r
357 #define FTOL(x,i,f) eifrac(x,i,f)
\r
359 /* i = &long integer input, x = &FLOAT output */
\r
360 #define LTOF(i,x) ltoe(i,x)
\r
362 /* Convert FLOAT a to decimal ASCII string with b digits */
\r
363 #define TOASC(a,b,c) etoasc(a,b,c)
\r
364 #endif /* not L128DOUBLE */
\r
368 /* The following subroutines are implementations of the above
\r
369 * named functions, using the native or default arithmetic.
\r
397 /* Important note: comparison can be done by subracting
\r
398 * or by a compare instruction that may or may not be
\r
399 * equivalent to subtracting.
\r
413 printf( "Compare fails\n" );
\r
447 *y = (FSIZE )ZFLOOR( *x );
\r
454 *y = (FSIZE )ZLOG( *x );
\r
461 *z = (FSIZE )ZPOW( *x, *y );
\r
468 *y = (FSIZE )ZSQRT( *x );
\r
479 y = (FSIZE )ZFLOOR( *x );
\r
509 sprintf( str, " %.17e ", x );
\r
512 /* default for FSETUP */
\r
516 #endif /* NATIVE */
\r
525 /*Small floating point constants.*/
\r
537 FLOAT(TwentySeven);
\r
543 /*Integer constants*/
\r
544 int NoTrials = 20; /*Number of tests for commutativity. */
\r
548 /* Definitions for declared types
\r
549 Guard == (Yes, No);
\r
550 Rounding == (Chopped, Rounded, Other);
\r
551 Message == packed array [1..40] of char;
\r
552 Class == (Flaw, Defect, Serious, Failure);
\r
564 typedef int Guard, Rounding, Class;
\r
565 typedef char Message;
\r
567 /* Declarations of Variables */
\r
627 int I, M, N, N1, stkflg;
\r
628 Guard GMult, GDiv, GAddSub;
\r
629 Rounding RMult, RDiv, RAddSub, RSqrt;
\r
630 int Break, Done, NotMonot, Monot, Anomaly, IEEE;
\r
631 int SqRWrng, UfNGrad;
\r
636 long lngint, lng2; /* intermediate for conversion between int and FLOAT */
\r
638 /* Computed constants. */
\r
639 /*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
\r
640 /*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
\r
649 /* Number of 16-bit groups to display */
\r
652 #define NPRT (sizeof( long double )/2)
\r
654 #define NPRT (sizeof( double )/2)
\r
661 printf( "%s\n", s );
\r
662 for( i=0; i<NPRT; i++ )
\r
663 printf( "%04x ", x[i] & 0xffff );
\r
667 /* define NOSIGNAL */
\r
669 #include <signal.h>
\r
671 #include <setjmp.h>
\r
673 /*typedef int (*Sig_type)();*/
\r
674 typedef void (*Sig_type)();
\r
677 /* Floating point exception receiver */
\r
681 printf( "\n* * * FLOATING-POINT ERROR * * *\n" );
\r
682 /* reinitialize the floating point unit */
\r
688 signal( SIGFPE, sigsave );
\r
691 longjmp( ovfl_buf, 1 );
\r
700 /* Do coprocessor or other initializations */
\r
704 "This version of paranoia omits test for extra precise subexpressions\n" );
\r
705 printf( "and includes a few additional tests.\n" );
\r
713 add( One, One, Two );
\r
714 printf( "1+1 = " );
\r
716 add( Two, One, Three );
\r
717 add( Three, One, Four );
\r
718 add( Four, One, Five );
\r
719 add( Five, One, Six );
\r
720 add( Four, Four, Eight );
\r
721 mul( Three, Three, Nine );
\r
722 add( Nine, One, Ten );
\r
723 mul( Nine, Three, TwentySeven );
\r
724 mul( Four, Eight, ThirtyTwo );
\r
725 mul( Four, Five, t );
\r
726 mul( t, Three, t );
\r
727 mul( t, Four, TwoForty );
\r
728 mov( One, MinusOne );
\r
730 div( Two, One, Half );
\r
731 add( One, Half, OneAndHalf );
\r
732 ErrCnt[Failure] = 0;
\r
733 ErrCnt[Serious] = 0;
\r
734 ErrCnt[Defect] = 0;
\r
738 signal( SIGFPE, sigfpe );
\r
740 printf("Program is now RUNNING tests on small integers:\n");
\r
742 add( Zero, Zero, t );
\r
743 if( cmp( t, Zero ) != 0)
\r
745 printf( "0+0 != 0\n" );
\r
746 ErrCnt[Failure] += 1;
\r
748 sub( One, One, t );
\r
749 if( cmp( t, Zero ) != 0 )
\r
751 printf( "1-1 != 0\n" );
\r
752 ErrCnt[Failure] += 1;
\r
754 if( cmp( One, Zero ) <= 0 )
\r
756 printf( "1 <= 0\n" );
\r
757 ErrCnt[Failure] += 1;
\r
759 add( One, One, t );
\r
760 if( cmp( t, Two ) != 0 )
\r
762 printf( "1+1 != 2\n" );
\r
763 ErrCnt[Failure] += 1;
\r
768 if( cmp(t,Zero) != 0 )
\r
770 ErrCnt[Serious] += 1;
\r
771 printf( "FLOOR(-0) should equal 0, is = " );
\r
774 if( cmp(Z, Zero) != 0)
\r
776 ErrCnt[Failure] += 1;
\r
777 printf("Comparison alleges that -0.0 is Non-zero!\n");
\r
781 div( TwoForty, One, U1 ); /* U1 = 0.001 */
\r
785 add( Two, One, t );
\r
786 if( cmp( t, Three ) != 0 )
\r
788 printf( "2+1 != 3\n" );
\r
789 ErrCnt[Failure] += 1;
\r
791 add( Three, One, t );
\r
792 if( cmp( t, Four ) != 0 )
\r
794 printf( "3+1 != 4\n" );
\r
795 ErrCnt[Failure] += 1;
\r
801 if( cmp( t, Zero ) != 0 )
\r
803 printf( "4+2*(-2) != 0\n" );
\r
804 ErrCnt[Failure] += 1;
\r
806 sub( Three, Four, t );
\r
808 if( cmp( t, Zero ) != 0 )
\r
810 printf( "4-3-1 != 0\n" );
\r
811 ErrCnt[Failure] += 1;
\r
813 sub( One, Zero, t );
\r
814 if( cmp( t, MinusOne ) != 0 )
\r
816 printf( "-1 != 0-1\n" );
\r
817 ErrCnt[Failure] += 1;
\r
819 add( One, MinusOne, t );
\r
820 if( cmp( t, Zero ) != 0 )
\r
822 printf( "1+(-1) != 0\n" );
\r
823 ErrCnt[Failure] += 1;
\r
827 add( MinusOne, t, t );
\r
828 if( cmp( t, Zero ) != 0 )
\r
830 printf( "-1+abs(1) != 0\n" );
\r
831 ErrCnt[Failure] += 1;
\r
833 mul( MinusOne, MinusOne, t );
\r
834 add( MinusOne, t, t );
\r
835 if( cmp( t, Zero ) != 0 )
\r
837 printf( "-1+(-1)*(-1) != 0\n" );
\r
838 ErrCnt[Failure] += 1;
\r
840 add( Half, MinusOne, t );
\r
842 if( cmp( t, Zero ) != 0 )
\r
844 printf( "1/2 + (-1) + 1/2 != 0\n" );
\r
845 ErrCnt[Failure] += 1;
\r
848 mul( Three, Three, t );
\r
849 if( cmp( t, Nine ) != 0 )
\r
851 printf( "3*3 != 9\n" );
\r
852 ErrCnt[Failure] += 1;
\r
854 mul( Nine, Three, t );
\r
855 if( cmp( t, TwentySeven ) != 0 )
\r
857 printf( "3*9 != 27\n" );
\r
858 ErrCnt[Failure] += 1;
\r
860 add( Four, Four, t );
\r
861 if( cmp( t, Eight ) != 0 )
\r
863 printf( "4+4 != 8\n" );
\r
864 ErrCnt[Failure] += 1;
\r
866 mul( Eight, Four, t );
\r
867 if( cmp( t, ThirtyTwo ) != 0 )
\r
869 printf( "8*4 != 32\n" );
\r
870 ErrCnt[Failure] += 1;
\r
872 sub( TwentySeven, ThirtyTwo, t );
\r
875 if( cmp( t, Zero ) != 0 )
\r
877 printf( "32-27-4-1 != 0\n" );
\r
878 ErrCnt[Failure] += 1;
\r
880 add( Four, One, t );
\r
881 if( cmp( t, Five ) != 0 )
\r
883 printf( "4+1 != 5\n" );
\r
884 ErrCnt[Failure] += 1;
\r
886 mul( Four, Five, t );
\r
887 mul( Three, t, t );
\r
889 if( cmp( t, TwoForty ) != 0 )
\r
891 printf( "4*5*3*4 != 240\n" );
\r
892 ErrCnt[Failure] += 1;
\r
894 div( Three, TwoForty, t );
\r
895 mul( Four, Four, t2 );
\r
896 mul( Five, t2, t2 );
\r
898 if( cmp( t, Zero ) != 0 )
\r
900 printf( "240/3 - 4*4*5 != 0\n" );
\r
901 ErrCnt[Failure] += 1;
\r
903 div( Four, TwoForty, t );
\r
904 mul( Five, Three, t2 );
\r
905 mul( Four, t2, t2 );
\r
907 if( cmp( t, Zero ) != 0 )
\r
909 printf( "240/4 - 5*3*4 != 0\n" );
\r
910 ErrCnt[Failure] += 1;
\r
912 div( Five, TwoForty, t );
\r
913 mul( Four, Three, t2 );
\r
914 mul( Four, t2, t2 );
\r
916 if( cmp( t, Zero ) != 0 )
\r
918 printf( "240/5 - 4*3*4 != 0\n" );
\r
919 ErrCnt[Failure] += 1;
\r
921 if(ErrCnt[Failure] == 0)
\r
923 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n\n");
\r
925 printf("Searching for Radix and Precision.\n");
\r
935 add( MinusOne, t, t );
\r
936 k = cmp( t, Zero );
\r
939 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
\r
940 mov( Zero, Precision );
\r
944 add( W, Y, Radix );
\r
946 sub( W, Radix, Radix );
\r
947 k = cmp( Radix, Zero );
\r
951 if( cmp(Radix, Two) < 0 )
\r
953 printf("Radix = " );
\r
955 if( cmp(Radix, One) != 0)
\r
960 add( One, Precision, Precision );
\r
961 mul( W, Radix, W );
\r
968 /* now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 */
\r
970 mul( Radix, U1, U2 );
\r
971 printf( "Closest relative separation found is U 1 = " );
\r
973 printf( "Recalculating radix and precision." );
\r
975 /*save old values*/
\r
979 mov( Precision, E3 );
\r
981 div( Three, Four, X );
\r
982 sub( One, X, Third );
\r
983 sub( Third, Half, F6 );
\r
985 sub( Third, X, X );
\r
987 if( cmp(X, U2) < 0 )
\r
990 /*... now X = (unknown no.) ulps of 1+...*/
\r
994 /* Y = Half * U2 + ThirtyTwo * U2 * U2; */
\r
995 mul( ThirtyTwo, U2, t );
\r
997 mul( Half, U2, Y );
\r
1002 k2 = cmp( X, Zero );
\r
1004 while ( ! ((k <= 0) || (k2 <= 0)));
\r
1006 /*... now U2 == 1 ulp of 1 + ... */
\r
1007 div( Three, Two, X );
\r
1008 sub( Half, X, F6 );
\r
1009 add( F6, F6, Third );
\r
1010 sub( Half, Third, X );
\r
1013 if( cmp(X, U1) < 0 )
\r
1016 /*... now X == (unknown no.) ulps of 1 -... */
\r
1020 /* Y = Half * U1 + ThirtyTwo * U1 * U1;*/
\r
1021 mul( ThirtyTwo, U1, t );
\r
1023 mul( Half, U1, Y );
\r
1025 sub( Y, Half, Y );
\r
1026 add( Half, Y, X );
\r
1027 sub( X, Half, Y );
\r
1028 add( Half, Y, X );
\r
1030 k2 = cmp( X, Zero );
\r
1031 } while ( ! ((k <= 0) || (k2 <= 0)));
\r
1032 /*... now U1 == 1 ulp of 1 - ... */
\r
1033 if( cmp( U1, E1 ) == 0 )
\r
1034 printf("confirms closest relative separation U1 .\n");
\r
1037 printf("gets better closest relative separation U1 = " );
\r
1040 div( U1, One, W );
\r
1041 sub( U1, Half, F9 );
\r
1042 add( F9, Half, F9 );
\r
1044 div( TwoForty, One, t2 );
\r
1046 FLOOR( t, Radix );
\r
1047 if( cmp(Radix, E0) == 0 )
\r
1048 printf("Radix confirmed.\n");
\r
1051 printf("MYSTERY: recalculated Radix = " );
\r
1055 add( Eight, Eight, t );
\r
1056 if( cmp( Radix, t ) > 0 )
\r
1058 printf( "Radix is too big: roundoff problems\n" );
\r
1059 ErrCnt[Defect] += 1;
\r
1062 if( cmp( Radix, Two ) == 0 )
\r
1064 if( cmp( Radix, Ten ) == 0 )
\r
1066 if( cmp( Radix, One ) == 0 )
\r
1070 printf( "Radix is not as good as 2 or 10\n" );
\r
1071 ErrCnt[Flaw] += 1;
\r
1073 /*=============================================*/
\r
1075 /*=============================================*/
\r
1076 sub( Half, F9, t );
\r
1077 if( cmp( t, Half ) >= 0 )
\r
1079 printf( "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?\n" );
\r
1080 ErrCnt[Failure] += 1;
\r
1084 sub( Half, X, Y );
\r
1085 sub( Half, Y, Z );
\r
1086 if( (cmp( X, One ) == 0) && (cmp( Z, Zero) != 0) )
\r
1088 printf( "Comparison is fuzzy ,X=1 but X-1/2-1/2 != 0\n" );
\r
1089 ErrCnt[Failure] += 1;
\r
1091 add( One, U2, X );
\r
1093 /*=============================================*/
\r
1095 /*=============================================*/
\r
1096 /*... BMinusU2 = nextafter(Radix, 0) */
\r
1098 sub( One, Radix, BMinusU2 );
\r
1099 sub( U2, BMinusU2, t );
\r
1100 add( One, t, BMinusU2 );
\r
1101 /* Purify Integers */
\r
1102 if( cmp(Radix,One) != 0 )
\r
1104 /*X = - TwoForty * LOG(U1) / LOG(Radix);*/
\r
1108 mul( TwoForty, X, X );
\r
1111 add( Half, X, Y );
\r
1115 mul( Four, t, t );
\r
1116 if( cmp( t, One ) < 0 )
\r
1118 div( TwoForty, X, Precision );
\r
1119 add( Half, Precision, Y );
\r
1121 sub( Y, Precision, t );
\r
1123 mul( TwoForty, t, t );
\r
1124 if( cmp( t, Half ) < 0 )
\r
1125 mov( Y, Precision );
\r
1127 FLOOR( Precision, t );
\r
1128 if( (cmp( Precision, t ) != 0) || (cmp( Radix, One ) == 0) )
\r
1130 printf("Precision cannot be characterized by an Integer number\n");
\r
1131 printf("of significant digits but, by itself, this is a minor flaw.\n");
\r
1133 if( cmp(Radix, One) == 0 )
\r
1134 printf("logarithmic encoding has precision characterized solely by U1.\n");
\r
1137 printf("The number of significant digits of the Radix is " );
\r
1138 show( Precision );
\r
1140 mul( U2, Nine, t );
\r
1141 mul( Nine, t, t );
\r
1142 mul( TwoForty, t, t );
\r
1143 if( cmp( t, One ) >= 0 )
\r
1145 printf( "Precision worse than 5 decimal figures\n" );
\r
1146 ErrCnt[Serious] += 1;
\r
1148 /*=============================================*/
\r
1150 /*=============================================*/
\r
1151 /* Test for extra-precise subepressions has been deleted. */
\r
1153 /*=============================================*/
\r
1154 if( cmp(Radix,Two) >= 0 )
\r
1156 mul( Radix, Radix, t );
\r
1162 if( cmp( X, U2 ) != 0 )
\r
1164 printf( "Subtraction is not normalized X=Y,X+Z != Y+Z!\n" );
\r
1165 ErrCnt[Failure] += 1;
\r
1167 if( cmp(X,U2) == 0 )
\r
1168 printf("Subtraction appears to be normalized, as it should be.");
\r
1171 printf("\nChecking for guard digit in *, /, and -.\n");
\r
1172 mul( F9, One, Y );
\r
1173 mul( One, F9, Z );
\r
1174 sub( Half, F9, X );
\r
1175 sub( Half, Y, Y );
\r
1177 sub( Half, Z, Z );
\r
1179 add( One, U2, X );
\r
1180 mul( X, Radix, T );
\r
1181 mul( Radix, X, R );
\r
1182 sub( Radix, T, X );
\r
1183 mul( Radix, U2, t );
\r
1185 sub( Radix, R, T );
\r
1186 mul( Radix, U2, t );
\r
1188 sub( One, Radix, t );
\r
1190 sub( One, Radix, t );
\r
1202 ErrCnt[Serious] += 1;
\r
1203 printf( "* lacks a Guard Digit, so 1*X != X\n" );
\r
1205 mul( Radix, U2, Z );
\r
1212 sub( U2, One, X );
\r
1218 if( (cmp(Y,Zero) > 0) || (cmp(Z,Zero) > 0) )
\r
1220 ErrCnt[Failure] += 1;
\r
1221 printf( "* gets too many final digits wrong.\n" );
\r
1223 sub( U2, One, Y );
\r
1224 add( One, U2, X );
\r
1227 div( Three, One, X );
\r
1228 div( Nine, Three, Z );
\r
1230 div( TwentySeven, Nine, T );
\r
1232 k = cmp( X, Zero );
\r
1233 k |= cmp( Y, Zero );
\r
1234 k |= cmp( Z, Zero );
\r
1237 ErrCnt[Defect] += 1;
\r
1238 printf( "Division lacks a Guard Digit, so error can exceed 1 ulp\n" );
\r
1239 printf( "or 1/3 and 3/9 and 9/27 may disagree\n" );
\r
1241 div( One, F9, Y );
\r
1242 sub( Half, F9, X );
\r
1243 sub( Half, Y, Y );
\r
1245 add( One, U2, X );
\r
1248 k = cmp( X, Zero );
\r
1249 k |= cmp( Y, Zero );
\r
1250 k |= cmp( Z, Zero );
\r
1256 ErrCnt[Serious] += 1;
\r
1257 printf( "Division lacks a Guard Digit, so X/1 != X\n" );
\r
1259 add( One, U2, X );
\r
1261 sub( Half, X, Y );
\r
1262 sub( Half, Y, Y );
\r
1263 if( cmp(Y,Zero) >= 0 )
\r
1265 ErrCnt[Serious] += 1;
\r
1266 printf( "Computed value of 1/1.000..1 >= 1\n" );
\r
1268 sub( U2, One, X );
\r
1269 mul( Radix, U2, Y );
\r
1271 mul( X, Radix, Z );
\r
1272 mul( Y, Radix, T );
\r
1273 div( Radix, Z, R );
\r
1274 div( Radix, T, StickyBit );
\r
1276 sub( Y, StickyBit, Y );
\r
1277 k = cmp( X, Zero );
\r
1278 k |= cmp( Y, Zero );
\r
1281 ErrCnt[Failure] += 1;
\r
1282 printf( "* and/or / gets too many last digits wrong\n" );
\r
1284 sub( U1, One, Y );
\r
1285 sub( F9, One, X );
\r
1287 sub( U2, Radix, T );
\r
1288 sub( BMinusU2, Radix, Z );
\r
1289 sub( T, Radix, T );
\r
1291 k |= cmp( Y, U1 );
\r
1292 k |= cmp( Z, U2 );
\r
1293 k |= cmp( T, U2 );
\r
1299 ErrCnt[Serious] += 1;
\r
1300 printf( "- lacks Guard Digit, so cancellation is obscured\n" );
\r
1302 sub( One, F9, t );
\r
1303 if( (cmp(F9,One) != 0) && (cmp(t,Zero) >= 0) )
\r
1305 ErrCnt[Serious] += 1;
\r
1306 printf("comparison alleges (1-U1) < 1 although\n");
\r
1307 printf(" subtration yields (1-U1) - 1 = 0 , thereby vitiating\n");
\r
1308 printf(" such precautions against division by zero as\n");
\r
1309 printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
\r
1311 if (GMult == Yes && GDiv == Yes && GAddSub == Yes)
\r
1312 printf(" *, /, and - appear to have guard digits, as they should.\n");
\r
1313 /*=============================================*/
\r
1315 /*=============================================*/
\r
1316 printf("Checking rounding on multiply, divide and add/subtract.\n");
\r
1320 div( Two, Radix, RadixD2 );
\r
1325 mov( Radix, AInvrse );
\r
1328 mov( AInvrse, X );
\r
1329 div( A1, AInvrse, AInvrse );
\r
1330 FLOOR( AInvrse, t );
\r
1331 k = cmp( t, AInvrse );
\r
1333 while( ! (k != 0 ) );
\r
1334 k = cmp( X, One );
\r
1335 k2 = cmp( A1, Three );
\r
1336 Done = (k == 0) || (k2 > 0);
\r
1338 add( Nine, One, A1 );
\r
1341 if( cmp(X, One) == 0 )
\r
1343 div( A1, One, AInvrse );
\r
1345 mov( AInvrse, Y );
\r
1350 sub( Half, Z, Z );
\r
1351 if( cmp( Z, Half ) != 0 )
\r
1353 ErrCnt[Failure] += 1;
\r
1354 printf( "X * (1/X) differs from 1\n" );
\r
1356 k = cmp( X, Radix );
\r
1363 add( One, U2, Y2 );
\r
1364 sub( U2, One, YY1 );
\r
1365 sub( U2, OneAndHalf, X );
\r
1366 add( OneAndHalf, U2, Y );
\r
1375 sub( OneAndHalf, X, X );
\r
1376 sub( OneAndHalf, Y, Y );
\r
1377 k = cmp( X, Zero );
\r
1378 k |= cmp( Y, Zero );
\r
1379 k |= cmp( Z, Zero );
\r
1380 if( cmp( T, Zero ) > 0 )
\r
1384 add( OneAndHalf, U2, X );
\r
1386 sub( U2, OneAndHalf, Y );
\r
1388 add( OneAndHalf, U2, Z );
\r
1390 sub( U2, OneAndHalf, T );
\r
1394 mul( Y, YY1, StickyBit );
\r
1398 add( StickyBit, Y, Y );
\r
1399 /* Z = S - (Z + U2 + U2); */
\r
1404 mul( t, YY1, StickyBit );
\r
1405 mul( Y2, YY1, YY1 );
\r
1406 sub( Y2, StickyBit, StickyBit );
\r
1407 sub( Half, YY1, YY1 );
\r
1408 k = cmp( X, Zero );
\r
1409 k |= cmp( Y, Zero );
\r
1410 k |= cmp( Z, Zero );
\r
1411 k |= cmp( T, Zero );
\r
1412 k |= cmp( StickyBit, Zero );
\r
1413 k |= cmp( YY1, Half );
\r
1417 printf("Multiplication appears to round correctly.\n");
\r
1422 k = cmp( t, Zero );
\r
1423 if( cmp( Y, Zero ) >= 0 )
\r
1426 k |= cmp( t, Zero );
\r
1427 if( cmp( T, Zero ) >= 0 )
\r
1429 add( StickyBit, U2, t );
\r
1430 k |= cmp( t, Zero );
\r
1431 if( cmp(YY1, Half) >= 0 )
\r
1435 printf("Multiplication appears to chop.\n");
\r
1439 printf("* is neither chopped nor correctly rounded.\n");
\r
1441 if( (RMult == Rounded) && (GMult == No) )
\r
1442 printf("Multiplication has inconsistent result");
\r
1446 printf("* is neither chopped nor correctly rounded.\n");
\r
1448 /*=============================================*/
\r
1450 /*=============================================*/
\r
1451 add( One, U2, Y2 );
\r
1452 sub( U2, One, YY1 );
\r
1453 add( OneAndHalf, U2, Z );
\r
1456 sub( U2, OneAndHalf, T );
\r
1462 sub( OneAndHalf, X, X );
\r
1465 add( OneAndHalf, U2, t );
\r
1467 sub( OneAndHalf, U2, t );
\r
1470 if( cmp( X, Zero ) > 0 )
\r
1472 if( cmp( Y, Zero ) > 0 )
\r
1474 if( cmp( Z, Zero ) > 0 )
\r
1476 if( cmp( T, Zero ) > 0 )
\r
1480 div( Y2, OneAndHalf, X );
\r
1481 sub( U2, OneAndHalf, Y );
\r
1482 add( U2, OneAndHalf, Z );
\r
1484 div( YY1, OneAndHalf, T );
\r
1490 add( Y2, U2, YY1 );
\r
1491 div( Y2, YY1, YY1 );
\r
1492 sub( OneAndHalf, Z, Z );
\r
1493 sub( Y2, YY1, Y2 );
\r
1494 sub( U1, F9, YY1 );
\r
1495 div( F9, YY1, YY1 );
\r
1496 k = cmp( X, Zero );
\r
1497 k |= cmp( Y, Zero );
\r
1498 k |= cmp( Z, Zero );
\r
1499 k |= cmp( T, Zero );
\r
1500 k |= cmp( Y2, Zero );
\r
1501 sub( Half, YY1, t );
\r
1502 sub( Half, F9, t2 );
\r
1503 k |= cmp( t, t2 );
\r
1507 printf("Division appears to round correctly.\n");
\r
1509 printf("Division test inconsistent\n");
\r
1514 if( cmp( X, Zero ) >= 0 )
\r
1516 if( cmp( Y, Zero ) >= 0 )
\r
1518 if( cmp( Z, Zero ) >= 0 )
\r
1520 if( cmp( T, Zero ) >= 0 )
\r
1522 if( cmp( Y, Zero ) >= 0 )
\r
1524 sub( Half, YY1, t );
\r
1525 sub( Half, F9, t2 );
\r
1526 if( cmp( t, t2 ) >= 0 )
\r
1531 printf("Division appears to chop.\n");
\r
1536 printf("/ is neither chopped nor correctly rounded.\n");
\r
1537 div( Radix, One, BInvrse );
\r
1538 mul( BInvrse, Radix, t );
\r
1539 sub( Half, t, t );
\r
1540 if( cmp( t, Half ) != 0 )
\r
1542 ErrCnt[Failure] += 1;
\r
1543 printf( "Radix * ( 1 / Radix ) differs from 1\n" );
\r
1547 /*=============================================*/
\r
1549 sub( Half, t, t );
\r
1550 k = cmp( t, Half );
\r
1551 add( BMinusU2, U2, t );
\r
1553 sub( One, Radix, t2 );
\r
1554 k |= cmp( t, t2 );
\r
1557 ErrCnt[Failure] += 1;
\r
1558 printf( "Incomplete carry-propagation in Addition\n" );
\r
1562 sub( U2, One, Y );
\r
1565 sub( Half, F9, Z );
\r
1566 sub( Half, X, X );
\r
1569 if( (cmp(X,Zero) == 0) && (cmp(Y,Zero) == 0) )
\r
1571 RAddSub = Chopped;
\r
1572 printf("Add/Subtract appears to be chopped.\n");
\r
1574 if(GAddSub == Yes)
\r
1576 add( Half, U2, X );
\r
1578 sub( U2, Half, Y );
\r
1582 add( One, U2, t );
\r
1587 printf( "1+U2-[u2(1/2+U2)+1] != 0\n" );
\r
1590 printf( "1-[U2(1/2-U2)+1] != 0\n" );
\r
1594 add( Half, U2, X );
\r
1596 sub( U2, Half, Y );
\r
1604 printf( "F9-[1-U1(1/2+U2)] != 0\n" );
\r
1607 printf( "1-[1-U1(1/2-U2)] != 0\n" );
\r
1611 RAddSub = Rounded;
\r
1612 printf("Addition/Subtraction appears to round correctly.\n");
\r
1614 printf( "Add/Subtract test inconsistent\n");
\r
1618 printf("Addition/Subtraction neither rounds nor chops.\n");
\r
1622 printf("Addition/Subtraction neither rounds nor chops.\n");
\r
1625 printf("Addition/Subtraction neither rounds nor chops.\n");
\r
1628 add( One, Half, X );
\r
1629 mul( Half, X, X );
\r
1631 add( One, U2, Y );
\r
1632 mul( Y, Half, Y );
\r
1635 add( Z, T, StickyBit );
\r
1636 if( cmp(StickyBit, Zero) != 0 )
\r
1639 ErrCnt[Flaw] += 1;
\r
1640 printf( "(X - Y) + (Y - X) is non zero!\n" );
\r
1642 mov( Zero, StickyBit );
\r
1643 FLOOR( RadixD2, t );
\r
1644 k2 = cmp( t, RadixD2 );
\r
1646 if( (GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
\r
1647 && (RMult == Rounded) && (RDiv == Rounded)
\r
1648 && (RAddSub == Rounded) && (k2 == 0) )
\r
1650 printf("Checking for sticky bit.\n");
\r
1652 add( Half, U1, X );
\r
1654 mul( Half, U2, Y );
\r
1658 sub( One, T, t2 );
\r
1659 if( cmp(t,Zero) > 0 )
\r
1662 printf( "[1+(1/2)U2]-1 > 0\n" );
\r
1664 if( cmp(t2,U2) < 0 )
\r
1667 printf( "[1+U2(1/2+U1)]-1 < U2\n" );
\r
1673 if( cmp(t,U2) < 0 )
\r
1676 printf( "[[1+U2(1/2+U1)]+(1/2)U2]-[1+U2(1/2+U1)] < U2\n" );
\r
1678 if( cmp(t2,Zero) != 0 )
\r
1681 printf( "(1/2)U2-[1+U2(1/2+U1)] != 0\n" );
\r
1683 add( Half, U1, X );
\r
1685 mul( Half, U1, Y );
\r
1690 if( cmp(t,Zero) != 0 )
\r
1693 printf( "(1-(1/2)U1)-1 != 0\n" );
\r
1695 if( cmp(t2,Zero) != 0 )
\r
1698 printf( "[1-U1(1/2+U1)]-F9 != 0\n" );
\r
1700 sub( U1, Half, Z );
\r
1705 if( cmp( t, Zero ) != 0 )
\r
1708 printf( "[F9-U1(1/2-U1)]-F9 != 0\n" );
\r
1712 if( cmp( t, Zero ) != 0 )
\r
1715 printf( "(F9-U1)-(F9-(1/2)U1) != 0\n" );
\r
1717 add( One, U2, Z );
\r
1718 mul( Z, OneAndHalf, Z );
\r
1719 add( OneAndHalf, U2, T );
\r
1722 div( Radix, Half, X );
\r
1724 mul( Radix, U2, Y );
\r
1727 if( cmp( T, Zero ) != 0 )
\r
1730 printf( "(3/2+U2)-3/2(1+U2)+U2 != 0\n" );
\r
1732 mul( Radix, U2, t );
\r
1735 if( cmp( t, Zero ) != 0 )
\r
1738 printf( "(1+1/2Radix)+Radix*U2-[1+1/(2Radix)][1+Radix*U2] != 0\n" );
\r
1740 if( cmp(Radix, Two) != 0 )
\r
1742 add( Two, U2, X );
\r
1745 if( cmp( t, Zero) != 0 )
\r
1751 printf("Sticky bit apparently used correctly.\n");
\r
1752 mov( One, StickyBit );
\r
1756 printf("Sticky bit used incorrectly or not at all.\n");
\r
1759 if( GMult == No || GDiv == No || GAddSub == No ||
\r
1760 RMult == Other || RDiv == Other || RAddSub == Other)
\r
1762 ErrCnt[Flaw] += 1;
\r
1763 printf("lack(s) of guard digits or failure(s) to correctly round or chop\n");
\r
1764 printf( "(noted above) count as one flaw in the final tally below\n" );
\r
1766 /*=============================================*/
\r
1768 /*=============================================*/
\r
1770 printf("Does Multiplication commute? ");
\r
1771 printf("Testing on %d random pairs.\n", NoTrials);
\r
1772 SQRT( Three, Random9 );
\r
1773 mov( Third, Random1 );
\r
1778 mov( Random1, X );
\r
1780 mov( Random1, Y );
\r
1786 while ( ! ((I > NoTrials) || (cmp(Z9,Zero) != 0)));
\r
1789 div( Three, Half, t );
\r
1790 add( One, t, Random1 );
\r
1792 add( t, One, Random2 );
\r
1793 mul( Random1, Random2, Z );
\r
1794 mul( Random2, Random1, Y );
\r
1795 /* Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
\r
1796 * Three) * ((U2 + U1) + One);
\r
1798 div( Three, Half, t2 );
\r
1799 add( One, t2, t2 );
\r
1806 if(! ((I == NoTrials) || (cmp(Z9,Zero) == 0)))
\r
1808 ErrCnt[Defect] += 1;
\r
1809 printf( "X * Y == Y * X trial fails.\n");
\r
1813 printf(" No failures found in %d integer pairs.\n", NoTrials);
\r
1815 /*=============================================*/
\r
1817 /*=============================================*/
\r
1824 printf("Seeking Underflow thresholds UfThold and E0.\n");
\r
1826 FLOOR( Precision, t );
\r
1827 if( cmp(Precision, t) != 0 )
\r
1829 mov( BInvrse, D );
\r
1830 mov( Precision, X );
\r
1833 mul( D, BInvrse, D );
\r
1836 while( cmp(X, Zero) > 0 );
\r
1840 /* ... D is power of 1/Radix < 1. */
\r
1842 if( setjmp(ovfl_buf) )
\r
1851 while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) );
\r
1859 if( setjmp(ovfl_buf) )
\r
1868 while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) );
\r
1873 if( cmp(Radix,Two) < 0 )
\r
1874 mov( Two, HInvrse );
\r
1876 mov( Radix, HInvrse );
\r
1877 div( HInvrse, One, H );
\r
1878 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
\r
1879 div( C, One, CInvrse );
\r
1882 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
\r
1884 if( setjmp(ovfl_buf) )
\r
1893 while( (cmp(E0,Z) > 0) && (cmp(t,Z) > 0) );
\r
1898 mov( E0, UfThold );
\r
1902 add( One, E9, S );
\r
1904 if( cmp(D,C) <= 0 )
\r
1906 mul( Radix, U2, E9 );
\r
1907 add( One, E9, S );
\r
1909 if( cmp(D, C) <= 0 )
\r
1911 ErrCnt[Failure] += 1;
\r
1912 printf( "multiplication gets too many last digits wrong.\n" );
\r
1913 mov( E0, Underflow );
\r
1915 mov( Z, PseudoZero );
\r
1920 mov( D, Underflow );
\r
1921 mul( Underflow, H, PseudoZero );
\r
1922 mov( Zero, UfThold );
\r
1925 mov( Underflow, YY1 );
\r
1926 mov( PseudoZero, Underflow );
\r
1928 if( cmp(t, E1) <= 0)
\r
1930 mul( Underflow, HInvrse, Y2 );
\r
1931 sub( Y2, YY1, E1 );
\r
1934 if( (cmp( UfThold, Zero ) == 0)
\r
1935 && (cmp(YY1, Y2) != 0) )
\r
1936 mov( YY1, UfThold );
\r
1938 mul( PseudoZero, H, PseudoZero );
\r
1939 add( PseudoZero, PseudoZero, t );
\r
1941 while( (cmp(Underflow, PseudoZero) > 0)
\r
1942 && (cmp(t, PseudoZero) > 0) );
\r
1944 /* Comment line 4530 .. 4560 */
\r
1945 if( cmp(PseudoZero, Zero) != 0 )
\r
1948 mov(PseudoZero, Z );
\r
1949 /* ... Test PseudoZero for "phoney- zero" violates */
\r
1950 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
\r
1952 if( cmp(PseudoZero, Zero) <= 0 )
\r
1954 ErrCnt[Failure] += 1;
\r
1955 printf("Positive expressions can underflow to an\n");
\r
1956 printf("allegedly negative value\n");
\r
1957 printf("PseudoZero that prints out as: " );
\r
1958 show( PseudoZero );
\r
1959 mov( PseudoZero, X );
\r
1961 if( cmp(X, Zero) <= 0 )
\r
1963 printf("But -PseudoZero, which should be\n");
\r
1964 printf("positive, isn't; it prints out as " );
\r
1970 ErrCnt[Flaw] += 1;
\r
1971 printf( "Underflow can stick at an allegedly positive\n");
\r
1972 printf("value PseudoZero that prints out as " );
\r
1973 show( PseudoZero );
\r
1978 /*=============================================*/
\r
1980 /*=============================================*/
\r
1981 mul( CInvrse, Y, t );
\r
1982 mul( CInvrse, YY1, t2 );
\r
1983 if( cmp(t,t2) > 0 )
\r
1986 mov( Underflow, E0 );
\r
1988 if(! ((cmp(E1,Zero) == 0) || (cmp(E1,E0) == 0)) )
\r
1990 ErrCnt[Defect] += 1;
\r
1991 if( cmp(E1,E0) < 0 )
\r
1993 printf("Products underflow at a higher");
\r
1994 printf(" threshold than differences.\n");
\r
1995 if( cmp(PseudoZero,Zero) == 0 )
\r
2000 printf("Difference underflows at a higher");
\r
2001 printf(" threshold than products.\n");
\r
2004 printf("Smallest strictly positive number found is E0 = " );
\r
2008 mov( E0, Underflow );
\r
2010 mov( Y, Underflow );
\r
2012 if( cmp(E1,Zero) == 0 )
\r
2014 if( cmp( UfThold,Zero) == 0 )
\r
2020 mov( Underflow, UfThold );
\r
2021 mul( CInvrse, Q, t );
\r
2022 mul( CInvrse, Y, t2 );
\r
2024 if( cmp( t, t2 ) != 0 )
\r
2026 mov( Y, UfThold );
\r
2027 ErrCnt[Failure] += 1;
\r
2028 printf( "Either accuracy deteriorates as numbers\n");
\r
2029 printf("approach a threshold = " );
\r
2031 printf(" coming down from " );
\r
2033 printf(" or else multiplication gets too many last digits wrong.\n");
\r
2038 ErrCnt[Failure] += 1;
\r
2039 printf( "Underflow confuses Comparison which alleges that\n");
\r
2040 printf("Q == Y while denying that |Q - Y| == 0; these values\n");
\r
2041 printf("print out as Q = " );
\r
2043 printf( ", Y = " );
\r
2047 printf ("|Q - Y| = " );
\r
2049 mov( Q, UfThold );
\r
2058 sub( t, UfThold, t );
\r
2060 if( (cmp(Q,UfThold) == 0) && (cmp(E1,E0) == 0)
\r
2061 && (cmp(t,E1) <= 0) )
\r
2064 printf("Underflow is gradual; it incurs Absolute Error =\n");
\r
2065 printf("(roundoff in UfThold) < E0.\n");
\r
2066 mul( E0, CInvrse, Y );
\r
2067 add( OneAndHalf, U2, t );
\r
2069 add( One, U2, X );
\r
2070 mul( CInvrse, X, X );
\r
2072 IEEE = (cmp(t,E0) == 0);
\r
2075 printf( "((CInvrse E0) (1.5+U2)) / (CInvrse (1+U2)) != E0\n" );
\r
2076 printf( "CInvrse = " );
\r
2078 printf( "E0 = " );
\r
2080 printf( "U2 = " );
\r
2086 printf( "Y/X = " );
\r
2094 div( UfThold, Underflow, R );
\r
2096 if( cmp(R,H) <= 0)
\r
2098 mul( R, UfThold, Z );
\r
2099 /* X = Z * (One + R * H * (One + H));*/
\r
2108 mov( UfThold, Z );
\r
2109 /*X = Z * (One + H * H * (One + H));*/
\r
2117 /* if(! ((cmp(X,Z) == 0) || (cmp(t,Zero) != 0)) )*/
\r
2118 if( (cmp(X,Z) != 0) && (cmp(t,Zero) == 0) )
\r
2120 /* ErrCnt[Flaw] += 1;*/
\r
2121 ErrCnt[Serious] += 1;
\r
2124 printf( "\tis not equal to Z = " );
\r
2126 /* sub( Z, X, Z9 );*/
\r
2127 printf("yet X - Z yields " );
\r
2129 printf("which compares equal to " );
\r
2131 printf(" Should this NOT signal Underflow, ");
\r
2132 printf("this is a SERIOUS DEFECT\nthat causes ");
\r
2133 printf("confusion when innocent statements like\n");;
\r
2134 printf(" if (X == Z) ... else");
\r
2135 printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
\r
2136 printf("encounter Division by Zero although actually\n");
\r
2137 printf("X / Z = 1 + " );
\r
2139 sub( Half, t, t );
\r
2140 sub( Half, t, t );
\r
2144 printf("The Underflow threshold is " );
\r
2146 printf( "below which calculation may suffer larger Relative error than" );
\r
2147 printf( " merely roundoff.\n");
\r
2148 mul( U1, U1, Y2 );
\r
2151 if( cmp( Y2,UfThold) <= 0 )
\r
2153 if( cmp(Y,E0) > 0 )
\r
2155 ErrCnt[Defect] += 1;
\r
2160 ErrCnt[Serious] += 1;
\r
2163 printf("Range is too narrow; U1^%d Underflows.\n", I);
\r
2167 /*Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;*/
\r
2168 LOG( UfThold, Y );
\r
2169 LOG( HInvrse, t );
\r
2171 mul( TwoForty, Y, Y );
\r
2172 sub( Y, Half, Y );
\r
2174 div( TwoForty, Y, Y );
\r
2176 sub( One, Y, Y2 ); /* ***** changed from Y2 = Y + Y */
\r
2177 printf("Since underflow occurs below the threshold\n");
\r
2178 printf("UfThold = " );
\r
2180 printf( "\tto the power " );
\r
2182 printf( "only underflow should afflict the expression " );
\r
2184 printf( "\tto the power " );
\r
2186 POW( HInvrse, Y2, V9 );
\r
2187 printf("Actually calculating yields: " );
\r
2189 add( Radix, Radix, t );
\r
2191 mul( t, UfThold, t );
\r
2192 if( (cmp(V9,Zero) < 0) || (cmp(V9,t) > 0) )
\r
2194 ErrCnt[Serious] += 1;
\r
2195 printf( "this is not between 0 and underflow\n");
\r
2196 printf(" threshold = " );
\r
2201 add( One, E9, t );
\r
2202 mul( UfThold, t, t );
\r
2203 if( cmp(V9,t) <= 0 )
\r
2204 printf("This computed value is O.K.\n");
\r
2207 ErrCnt[Defect] += 1;
\r
2208 printf( "this is not between 0 and underflow\n");
\r
2209 printf(" threshold = " );
\r
2218 /*=============================================*/
\r
2220 /*=============================================*/
\r
2222 printf("Searching for Overflow threshold:\n");
\r
2223 printf("This may generate an error.\n");
\r
2226 mov( CInvrse, Y ); /* a large power of 2 */
\r
2228 mul( HInvrse, Y, V9 ); /* HInvrse = 2 */
\r
2229 if (setjmp(ovfl_buf))
\r
2235 mul( HInvrse, Y, V9 );
\r
2237 while( cmp(V9,Y) < 0 ); /* V9 = 2 * Y */
\r
2243 printf( "\ttimes " );
\r
2245 printf( "\tequals " );
\r
2249 printf("Can `Z = -Y' overflow?\n");
\r
2250 printf("Trying it on Y = " );
\r
2257 if( cmp(t,t2) == 0 )
\r
2258 printf("Seems O.K.\n");
\r
2261 printf("finds a Flaw, -(-Y) differs from Y.\n");
\r
2262 printf( "V-Y=t:" );
\r
2266 printf( "V+V0=t2:" );
\r
2270 ErrCnt[Flaw] += 1;
\r
2272 if( (cmp(Z, Y) != 0) && (I != 0) )
\r
2274 ErrCnt[Serious] += 1;
\r
2275 printf("overflow past " );
\r
2277 printf( "\tshrinks to " );
\r
2279 printf( "= Y * " );
\r
2282 /*Y = V * (HInvrse * U2 - HInvrse);*/
\r
2283 mul( HInvrse, U2, Y );
\r
2284 sub( HInvrse, Y, Y );
\r
2286 /*Z = Y + ((One - HInvrse) * U2) * V;*/
\r
2287 sub( HInvrse, One, Z );
\r
2291 if( cmp(Z,V0) < 0 )
\r
2293 if( cmp(Y,V0) < 0)
\r
2296 if( cmp(t,V0) < 0 )
\r
2298 printf("Overflow threshold is V = " );
\r
2302 printf("Overflow saturates at V0 = " );
\r
2306 printf("There is no saturation value because the system traps on overflow.\n");
\r
2308 mul( V, One, V9 );
\r
2309 printf("No Overflow should be signaled for V * 1 = " );
\r
2311 div( One, V, V9 );
\r
2312 printf(" nor for V / 1 = " );
\r
2314 printf("Any overflow signal separating this * from the one\n");
\r
2315 printf("above is a DEFECT.\n");
\r
2316 /*=============================================*/
\r
2318 /*=============================================*/
\r
2322 if( cmp(t,V) >= 0 )
\r
2326 if( cmp(t,V0) >= 0 )
\r
2328 mov( UfThold, t );
\r
2330 if( cmp(t,V) >= 0 )
\r
2332 if( cmp(UfThold,V) >= 0 )
\r
2336 ErrCnt[Failure] += 1;
\r
2337 printf( "Comparisons involving +-");
\r
2341 printf("are confused by Overflow." );
\r
2343 /*=============================================*/
\r
2345 /*=============================================*/
\r
2347 for(Indx = 1; Indx <= 3; ++Indx) {
\r
2350 case 1: mov(UfThold, Z); break;
\r
2351 case 2: mov( E0, Z); break;
\r
2352 case 3: mov(PseudoZero, Z); break;
\r
2354 if( cmp(Z, Zero) != 0 )
\r
2358 mul( Radix, E9, t );
\r
2361 add( One, Radix, t2 );
\r
2362 add( t2, E9, t2 );
\r
2364 if( (cmp(t,Z) < 0) || (cmp(Y,t2) > 0) )
\r
2366 if( cmp(V9,U1) > 0 )
\r
2367 ErrCnt[Serious] += 1;
\r
2369 ErrCnt[Defect] += 1;
\r
2370 printf("Comparison alleges that what prints as Z = " );
\r
2372 printf(" is too far from sqrt(Z) ^ 2 = " );
\r
2380 for(Indx = 1; Indx <= 2; ++Indx)
\r
2387 mul( Radix, E9, X );
\r
2391 mul( Two, Radix, t );
\r
2395 if( (cmp(V9,t) < 0) || (cmp(V9,Z) > 0) )
\r
2398 if( cmp(X,W) < 0 )
\r
2399 ErrCnt[Serious] += 1;
\r
2401 ErrCnt[Defect] += 1;
\r
2402 printf("Comparison alleges that Z = " );
\r
2404 printf(" is too far from sqrt(Z) ^ 2 :" );
\r
2412 mul( UfThold, V, X );
\r
2413 mul( Radix, Radix, Y );
\r
2415 if( (cmp(t,One) < 0) || (cmp(X,Y) > 0) )
\r
2419 if( (cmp(t,U1) < 0) || (cmp(X,t2) > 0) )
\r
2421 ErrCnt[Defect] += 1;
\r
2422 printf( "Badly " );
\r
2426 ErrCnt[Flaw] += 1;
\r
2428 printf(" unbalanced range; UfThold * V = " );
\r
2430 printf( "\tis too far from 1.\n");
\r
2434 for(Indx = 1; Indx <= 5; ++Indx)
\r
2439 case 2: add( One, U2, X ); break;
\r
2440 case 3: mov( V, X ); break;
\r
2441 case 4: mov(UfThold,X); break;
\r
2442 case 5: mov(Radix,X);
\r
2447 if (setjmp(ovfl_buf))
\r
2449 printf(" X / X traps when X = " );
\r
2454 /*V9 = (Y / X - Half) - Half;*/
\r
2456 sub( Half, t, t );
\r
2457 sub( Half, t, V9 );
\r
2458 if( cmp(V9,Zero) == 0 )
\r
2462 if( (cmp(V9,t) == 0) && (Indx < 5) )
\r
2464 ErrCnt[Flaw] += 1;
\r
2468 ErrCnt[Serious] += 1;
\r
2470 printf(" X / X differs from 1 when X = " );
\r
2472 printf(" instead, X / X - 1/2 - 1/2 = " );
\r
2480 static char *msg[] = {
\r
2481 "FAILUREs encountered =",
\r
2482 "SERIOUS DEFECTs discovered =",
\r
2483 "DEFECTs discovered =",
\r
2484 "FLAWs discovered =" };
\r
2486 for(i = 0; i < 4; i++) if (ErrCnt[i])
\r
2487 printf("The number of %-29s %d.\n",
\r
2488 msg[i], ErrCnt[i]);
\r
2491 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
\r
2492 + ErrCnt[Flaw]) > 0) {
\r
2493 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
\r
2494 Defect] == 0) && (ErrCnt[Flaw] > 0)) {
\r
2495 printf("The arithmetic diagnosed seems ");
\r
2496 printf("satisfactory though flawed.\n");
\r
2498 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
\r
2499 && ( ErrCnt[Defect] > 0)) {
\r
2500 printf("The arithmetic diagnosed may be acceptable\n");
\r
2501 printf("despite inconvenient Defects.\n");
\r
2503 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
\r
2504 printf("The arithmetic diagnosed has ");
\r
2505 printf("unacceptable serious defects.\n");
\r
2507 if (ErrCnt[Failure] > 0) {
\r
2508 printf("Fatal FAILURE may have spoiled this");
\r
2509 printf(" program's subsequent diagnoses.\n");
\r
2513 printf("No failures, defects nor flaws have been discovered.\n");
\r
2514 if (! ((RMult == Rounded) && (RDiv == Rounded)
\r
2515 && (RAddSub == Rounded) && (RSqrt == Rounded)))
\r
2516 printf("The arithmetic diagnosed seems satisfactory.\n");
\r
2519 if( cmp( Radix, Two ) == 0 )
\r
2521 if( cmp( Radix, Ten ) == 0 )
\r
2523 if( (cmp(StickyBit,One) >= 0) && (k == 1) )
\r
2525 printf("Rounding appears to conform to ");
\r
2526 printf("the proposed IEEE standard P");
\r
2528 k |= cmp( Radix, Two );
\r
2529 mul( Four, Three, t );
\r
2531 sub( t, Precision, t );
\r
2532 sub( TwentySeven, Precision, t2 );
\r
2533 sub( TwentySeven, t2, t2 );
\r
2534 add( t2, One, t2 );
\r
2536 if( (cmp(Radix,Two) == 0)
\r
2537 && (cmp(t,Zero) == 0) )
\r
2545 printf(",\nexcept for possibly Double Rounding");
\r
2546 printf(" during Gradual Underflow.\n");
\r
2549 printf("The arithmetic diagnosed appears to be excellent!\n");
\r
2553 printf("\nA total of %d floating point exceptions were registered.\n",
\r
2555 printf("END OF TEST.\n");
\r
2560 /* Random computes
\r
2561 X = (Random1 + Random9)^5
\r
2562 Random1 = X - FLOOR(X) + 0.000005 * X;
\r
2563 and returns the new value of Random1
\r
2567 static int randflg = 0;
\r
2573 if( randflg == 0 )
\r
2578 mul( Five, t, C5em6 );
\r
2581 add( Random1, Random9, t );
\r
2583 mul( t2, t2, t2 );
\r
2587 mul( t, C5em6, t );
\r
2588 add( t, t2, Random1 );
\r
2589 /*return(Random1);*/
\r
2594 SqXMinX( ErrKind )
\r
2597 mul( X, BInvrse, t2 );
\r
2599 /*SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;*/
\r
2600 mul( X, X, Sqarg );
\r
2601 SQRT( Sqarg, SqEr );
\r
2602 sub( t2, SqEr, SqEr );
\r
2603 sub( t, SqEr, SqEr );
\r
2604 div( OneUlp, SqEr, SqEr );
\r
2605 if( cmp(SqEr,Zero) != 0)
\r
2609 ErrCnt[ErrKind] += 1;
\r
2610 printf("sqrt of " );
\r
2613 printf( "minus " );
\r
2615 printf( "equals " );
\r
2616 mul( OneUlp, SqEr, t );
\r
2618 printf("\tinstead of correct value 0 .\n");
\r
2627 /*X = FLOOR(Half - X / Radix) * Radix + X;*/
\r
2628 div( Radix, X, t );
\r
2629 sub( t, Half, t );
\r
2631 mul( t, Radix, t );
\r
2633 /*Q = (Q - X * Z) / Radix + X * X * (D / Radix);*/
\r
2636 div( Radix, t, t );
\r
2637 div( Radix, D, t2 );
\r
2641 /*Z = Z - Two * X * D;*/
\r
2646 if( cmp(Z,Zero) <= 0)
\r
2651 mul( Radix, D, D );
\r
2658 sub( Radix, X, t );
\r
2659 sub( Radix, Z2, t2 );
\r
2661 if( cmp(t,t2) < 0 )
\r
2665 if( cmp(t,t2) > 0 )
\r
2667 /*if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {*/
\r
2674 /*Y2 = (X2 - Z2) - (Y - Z2);*/
\r
2675 sub( Z2, X2, Y2 );
\r
2678 sub( Half, Y, X2 );
\r
2679 div( X2, X8, X2 );
\r
2680 mul( Half, X2, t );
\r
2683 /*SqEr = (Y2 + Half) + (Half - X2);*/
\r
2684 add( Y2, Half, SqEr );
\r
2685 sub( X2, Half, t );
\r
2686 add( t, SqEr, SqEr );
\r
2688 sub( X2, Y2, SqEr );
\r
2697 if( cmp(Y,X) != 0 )
\r
2701 if( (cmp(Z,Zero) == 0) && (cmp(Q,Zero) <= 0) )
\r
2702 printf("WARNING: computing\n");
\r
2705 ErrCnt[Defect] += 1;
\r
2706 printf( "computing\n");
\r
2709 printf( "\tto the power " );
\r
2711 printf("\tyielded " );
\r
2713 printf("\twhich compared unequal to correct " );
\r
2716 printf("\t\tthey differ by " );
\r
2719 N = N + 1; /* ... count discrepancies. */
\r
2731 /*Q = (FLOAT) I;*/
\r
2740 while( cmp(X,W) < 0 );
\r
2743 /* PrintIfNPositive */
\r
2745 PrintIfNPositive()
\r
2748 printf("Similar discrepancies have occurred %d times.\n", N);
\r
2757 if( cmp(Z,Zero) != 0)
\r
2761 printf("Since comparison denies Z = 0, evaluating ");
\r
2762 printf("(Z + Z) / Z should be safe.\n");
\r
2764 if (setjmp(ovfl_buf))
\r
2765 goto very_serious;
\r
2768 printf("What the machine gets for (Z + Z) / Z is " );
\r
2770 sub( Two, Q9, t );
\r
2772 mul( Radix, U2, t2 );
\r
2773 if( cmp(t,t2) < 0 )
\r
2775 printf("This is O.K., provided Over/Underflow");
\r
2776 printf(" has NOT just been signaled.\n");
\r
2780 if( (cmp(Q9,One) < 0) || (cmp(Q9,Two) > 0) )
\r
2784 ErrCnt [Serious] = ErrCnt [Serious] + 1;
\r
2785 printf("This is a VERY SERIOUS DEFECT!\n");
\r
2790 ErrCnt[Defect] += 1;
\r
2791 printf("This is a DEFECT!\n");
\r
2794 mul( Z, One, V9 );
\r
2795 mov( V9, Random1 );
\r
2796 mul( One, Z, V9 );
\r
2797 mov( V9, Random2 );
\r
2798 div( One, Z, V9 );
\r
2799 if( (cmp(Z,Random1) == 0) && (cmp(Z,Random2) == 0)
\r
2800 && (cmp(Z,V9) == 0) )
\r
2808 ErrCnt[Defect] += 1;
\r
2809 printf( "What prints as Z = ");
\r
2811 printf( "\tcompares different from " );
\r
2812 if( cmp(Z,Random1) != 0)
\r
2814 printf("Z * 1 = " );
\r
2817 if( (cmp(Z,Random2) != 0)
\r
2818 || (cmp(Random2,Random1) != 0) )
\r
2820 printf("1 * Z == " );
\r
2823 if( cmp(Z,V9) != 0 )
\r
2825 printf("Z / 1 = " );
\r
2828 if( cmp(Random2,Random1) != 0 )
\r
2830 ErrCnt[Defect] += 1;
\r
2831 printf( "Multiplication does not commute!\n");
\r
2832 printf("\tComparison alleges that 1 * Z = " );
\r
2834 printf("\tdiffers from Z * 1 = " );
\r
2850 if( cmp( x, Zero ) < 0 )
\r
2863 printf("\nRunning test of square root(x).\n");
\r
2868 k |= cmp( Zero, t );
\r
2872 k |= cmp( t, t2 );
\r
2874 k |= cmp( One, t );
\r
2877 ErrCnt[Failure] += 1;
\r
2878 printf( "Square root of 0.0, -0.0 or 1.0 wrong\n");
\r
2880 mov( Zero, MinSqEr );
\r
2881 mov( Zero, MaxSqEr );
\r
2884 mov( U2, OneUlp );
\r
2885 SqXMinX( Serious );
\r
2886 mov( BInvrse, X );
\r
2887 mul( BInvrse, U1, OneUlp );
\r
2888 SqXMinX( Serious );
\r
2890 mul( U1, U1, OneUlp );
\r
2891 SqXMinX( Serious );
\r
2892 if( cmp(J,Zero) != 0)
\r
2894 printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
\r
2898 if( cmp(Radix,One) != 0 )
\r
2900 lngint = NoTrials;
\r
2901 LTOF( &lngint, t );
\r
2902 FTOL( t, &lng2, X );
\r
2903 if( lngint != lng2 )
\r
2905 printf( "Integer conversion error\n" );
\r
2911 mul( Radix, Y, Y );
\r
2914 while( ! (cmp(t2,t) >= 0) );
\r
2916 mul( X, U2, OneUlp );
\r
2921 SqXMinX( Defect );
\r
2922 if( cmp(J,Zero) > 0 )
\r
2926 printf("Test for sqrt monotonicity.\n");
\r
2928 mov( BMinusU2, X );
\r
2930 mul( Radix, U2, Z );
\r
2931 add( Radix, Z, Z );
\r
2934 while( ! (NotMonot || Monot))
\r
2940 if( (cmp(X,Q) > 0) || (cmp(Q,Z) > 0) )
\r
2944 add( Q, Half, Q );
\r
2947 if( (I > 0) || (cmp(Radix,t) == 0) )
\r
2955 mul( Y, BInvrse, Y );
\r
2969 printf("sqrt has passed a test for Monotonicity.\n");
\r
2972 ErrCnt[Defect] += 1;
\r
2973 printf("sqrt(X) is non-monotonic for X near " );
\r
2976 /*=============================================*/
\r
2978 /*=============================================*/
\r
2979 add( MinSqEr, Half, MinSqEr );
\r
2980 sub( Half, MaxSqEr, MaxSqEr);
\r
2981 /*Y = (SQRT(One + U2) - One) / U2;*/
\r
2982 add( One, U2, Sqarg );
\r
2986 /*SqEr = (Y - One) + U2 / Eight;*/
\r
2988 div( Eight, U2, SqEr );
\r
2989 add( t, SqEr, SqEr );
\r
2991 div( Eight, U2, SqEr );
\r
2992 add( Y, SqEr, SqEr );
\r
2994 /*Y = ((SQRT(F9) - U2) - (One - U2)) / U1;*/
\r
2998 sub( U2, One, t );
\r
3001 div( Eight, U1, SqEr );
\r
3002 add( Y, SqEr, SqEr );
\r
3004 /*SqEr = (Y + One) + U1 / Eight;*/
\r
3005 div( Eight, U1, t );
\r
3006 add( Y, One, SqEr );
\r
3007 add( SqEr, t, SqEr );
\r
3009 mov( U2, OneUlp );
\r
3011 for( Indx = 1; Indx <= 3; ++Indx)
\r
3013 /*Y = SQRT((X + U1 + X) + F9);*/
\r
3019 /*Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;*/
\r
3020 sub( U2, One, t );
\r
3024 div( OneUlp, Y, Y );
\r
3025 /*Z = ((U1 - X) + F9) * Half * X * X / OneUlp;*/
\r
3028 mul( t, Half, t );
\r
3031 div( OneUlp, t, Z );
\r
3032 add( Y, Half, SqEr );
\r
3033 add( SqEr, Z, SqEr );
\r
3035 sub( Half, Y, SqEr );
\r
3036 add( SqEr, Z, SqEr );
\r
3038 if(((Indx == 1) || (Indx == 3)))
\r
3040 /*X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));*/
\r
3041 mov( OneUlp, Sqarg );
\r
3043 mul( Nine, t, t );
\r
3044 div( t, Eight, t );
\r
3048 mul( OneUlp, t, X );
\r
3052 mov( U1, OneUlp );
\r
3057 /*=============================================*/
\r
3059 /*=============================================*/
\r
3062 if( cmp(Radix,One) != 0 )
\r
3064 printf("Testing whether sqrt is rounded or chopped.\n");
\r
3065 /*D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));*/
\r
3066 FLOOR( Precision, t2 );
\r
3067 add( One, Precision, t );
\r
3069 POW( Radix, t, D );
\r
3070 add( Half, D, D );
\r
3072 /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
\r
3073 div( Radix, D, X );
\r
3077 if( (cmp(X,t) != 0) || (cmp(Y,t2) != 0) )
\r
3080 printf( "Anomaly 1\n" );
\r
3088 sub( One, Radix, Z1 );
\r
3089 mul( Four, D, FourD );
\r
3092 if( cmp(Y2,Z2) >0 )
\r
3098 /*X1 = FABS(Q + FLOOR(Half - Q / YY1) * YY1);*/
\r
3100 sub( t, Half, t );
\r
3108 while( ! (cmp(X1,Zero) <= 0) );
\r
3109 if( cmp(Q,One) <= 0 )
\r
3116 add( X, Eight, X );
\r
3118 if( cmp(Y2,FourD) >= 0 )
\r
3119 sub( FourD, Y2, Y2 );
\r
3121 while( ! (cmp(Y,D) >= 0) );
\r
3122 sub( Z2, FourD, X8 );
\r
3125 div( FourD, Q, Q );
\r
3126 div( Eight, X8, X8 );
\r
3128 if( cmp(Q,t) != 0 )
\r
3131 printf( "Anomaly 2\n" );
\r
3139 /*X = X - FLOOR(X / Radix) * Radix;*/
\r
3140 div( Radix, X, t );
\r
3142 mul( t, Radix, t );
\r
3144 if( cmp(X,One) == 0 )
\r
3147 sub( One, Z1, Z1 );
\r
3149 while( ! (Break || (cmp(Z1,Zero) <= 0)) );
\r
3150 if( (cmp(Z1,Zero) <= 0) && (! Break))
\r
3152 printf( "Anomaly 3\n" );
\r
3157 if( cmp(Z1,RadixD2) > 0)
\r
3158 sub( Radix, Z1, Z1 );
\r
3164 while( ! (cmp(t,F9) >= 0) );
\r
3165 mul( D, Radix, t );
\r
3168 if (cmp(t,t2) != 0 )
\r
3170 printf( "Anomaly 4\n" );
\r
3178 mul( t, Half, t );
\r
3184 mul( t, Half, t );
\r
3195 if(cmp(t,t2) != 0 )
\r
3197 printf( "Anomaly 5\n" );
\r
3202 /*Y = (D - Z2) + (Z2 + (One - Z) * Half);*/
\r
3204 mul( t, Half, t );
\r
3208 /*X = (D - Z2) + (Z2 - Z + Q);*/
\r
3215 mul( Y, Half, Y );
\r
3220 printf( "Anomaly 6\n" );
\r
3228 if ((I == 0) || Anomaly)
\r
3230 ErrCnt[Failure] += 1;
\r
3231 printf( "Anomalous arithmetic with Integer < \n");
\r
3232 printf("Radix^Precision = " );
\r
3234 printf(" fails test whether sqrt rounds or chops.\n");
\r
3240 if(! ((cmp(MinSqEr,Zero) < 0) || (cmp(MaxSqEr,Zero) > 0))) {
\r
3242 printf("Square root appears to be correctly rounded.\n");
\r
3247 add( MaxSqEr, U2, t );
\r
3248 sub( Half, U2, t2 );
\r
3249 if( cmp(t,t2) > 0 )
\r
3251 if( cmp( MinSqEr, Half ) > 0 )
\r
3253 add( MinSqEr, Radix, t );
\r
3254 if( cmp( t, Half ) < 0 )
\r
3261 printf("Square root appears to be chopped.\n");
\r
3267 printf("Square root is neither chopped nor correctly rounded.\n");
\r
3268 printf("Observed errors run from " );
\r
3269 sub( Half, MinSqEr, t );
\r
3272 add( Half, MaxSqEr, t );
\r
3274 printf( "ulps.\n" );
\r
3275 sub( MinSqEr, MaxSqEr, t );
\r
3276 mul( Radix, Radix, t2 );
\r
3277 if( cmp( t, t2 ) >= 0 )
\r
3279 ErrCnt[Serious] += 1;
\r
3280 printf( "sqrt gets too many last digits wrong\n");
\r
3292 if( cmp(SqEr,MinSqEr) < 0 )
\r
3295 mov( SqEr, MinSqEr );
\r
3300 if( cmp(SqEr,MaxSqEr) > 0 )
\r
3303 mov( SqEr, MaxSqEr );
\r
3309 printf( "Square root of " );
\r
3311 printf( "\tis in error by " );
\r
3321 /*=============================================*/
\r
3323 /*=============================================*/
\r
3325 printf("Testing powers Z^i for small Integers Z and i.\n");
\r
3327 /* ... test powers of zero. */
\r
3342 if( cmp(Z,MinusOne) == 0 )
\r
3346 mov( MinusOne, Z );
\r
3347 PrintIfNPositive();
\r
3349 /* .. if(-1)^N is invalid, replace MinusOne by One. */
\r
3354 PrintIfNPositive();
\r
3358 /*M = FLOOR(Two * LOG(W) / LOG(A1));*/
\r
3364 FTOL( t, &lngint, t2 );
\r
3372 if( cmp(Z,AInvrse) == 0 )
\r
3375 mov( AInvrse, Z );
\r
3377 while( ! (Break) );
\r
3378 /*=============================================*/
\r
3380 /*=============================================*/
\r
3381 /* Powers of Radix have been tested, */
\r
3382 /* next try a few primes */
\r
3395 div( Three, Z, t );
\r
3397 mul( Three, t, t );
\r
3399 while( cmp(t,Z) == 0 );
\r
3400 mul( Eight, Three, t );
\r
3402 while( cmp(Z,t) < 0 );
\r
3406 printf("Errors like this may invalidate financial calculations\n");
\r
3407 printf("\tinvolving interest rates.\n");
\r
3409 PrintIfNPositive();
\r
3412 printf("... no discrepancies found.\n");
\r
3415 else printf("\n");
\r
3423 /* ...calculate Exp2 == exp(2) == 7.38905 60989 30650 22723 04275-... */
\r
3425 mov( Two, t2 ); /*I = 2;*/
\r
3427 mul( Two, Three, Y );
\r
3433 add( t2, One, t2 ); /*I = I + 1;*/
\r
3435 div( t, Y, Y ); /*Y = Y / (I + I);*/
\r
3441 while( cmp(X,Z) > 0 );
\r
3443 /*Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);*/
\r
3444 div( Eight, One, t );
\r
3445 add( OneAndHalf, t, Z );
\r
3446 mul( OneAndHalf, ThirtyTwo, t );
\r
3450 mul( X, X, Exp2 );
\r
3453 printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = " );
\r
3455 printf( "\tas X -> 1.\n" );
\r
3458 sub( BInvrse, X, Z );
\r
3459 /*Z = (X + One) / (Z - (One - BInvrse));*/
\r
3460 add( X, One, t2 );
\r
3461 sub( BInvrse, One, t );
\r
3464 POW( X, Z, Sqarg );
\r
3465 sub( Exp2, Sqarg, Q );
\r
3468 mul( TwoForty, U2, t2 );
\r
3469 if( cmp( t, t2 ) > 0 )
\r
3472 sub( BInvrse, X, V9 );
\r
3473 sub( BInvrse, One, t );
\r
3475 ErrCnt[Defect] += 1;
\r
3476 printf( "Calculated " );
\r
3478 printf(" for \t(1 + " );
\r
3480 printf( "\tto the power " );
\r
3482 printf("\tdiffers from correct value by " );
\r
3484 printf("\tThis much error may spoil financial\n");
\r
3485 printf("\tcalculations involving tiny interest rates.\n");
\r
3498 if( (cmp(Z,One) > 0) && (I < NoTrials) )
\r
3502 if( cmp(X,One) > 0 )
\r
3505 printf("Accuracy seems adequate.\n");
\r
3510 add( One, U2, X );
\r
3518 /*=============================================*/
\r
3520 /*=============================================*/
\r
3521 printf("Testing powers Z^Q at four nearly extreme values.\n");
\r
3524 /*Q = FLOOR(Half - LOG(C) / LOG(A1));*/
\r
3528 sub( t, Half, t );
\r
3533 mov( CInvrse, X );
\r
3540 if( cmp(Z,One) < 0 )
\r
3543 mov( AInvrse, Z );
\r
3545 while( ! (Break));
\r
3546 PrintIfNPositive();
\r
3548 printf(" ... no discrepancies found.\n");
\r