URL
https://opencores.org/ocsvn/openrisc/openrisc/trunk
Subversion Repositories openrisc
[/] [openrisc/] [trunk/] [gnu-old/] [gcc-4.2.2/] [contrib/] [paranoia.cc] - Rev 868
Go to most recent revision | Compare with Previous | Blame | View Log
/* A C version of Kahan's Floating Point Test "Paranoia" Thos Sumner, UCSF, Feb. 1985 David Gay, BTL, Jan. 1986 This is a rewrite from the Pascal version by B. A. Wichmann, 18 Jan. 1985 (and does NOT exhibit good C programming style). Adjusted to use Standard C headers 19 Jan. 1992 (dmg); (C) Apr 19 1983 in BASIC version by: Professor W. M. Kahan, 567 Evans Hall Electrical Engineering & Computer Science Dept. University of California Berkeley, California 94720 USA converted to Pascal by: B. A. Wichmann National Physical Laboratory Teddington Middx TW11 OLW UK converted to C by: David M. Gay and Thos Sumner AT&T Bell Labs Computer Center, Rm. U-76 600 Mountain Avenue University of California Murray Hill, NJ 07974 San Francisco, CA 94143 USA USA with simultaneous corrections to the Pascal source (reflected in the Pascal source available over netlib). [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.] Reports of results on various systems from all the versions of Paranoia are being collected by Richard Karpinski at the same address as Thos Sumner. This includes sample outputs, bug reports, and criticisms. You may copy this program freely if you acknowledge its source. Comments on the Pascal version to NPL, please. The following is from the introductory commentary from Wichmann's work: The BASIC program of Kahan is written in Microsoft BASIC using many facilities which have no exact analogy in Pascal. The Pascal version below cannot therefore be exactly the same. Rather than be a minimal transcription of the BASIC program, the Pascal coding follows the conventional style of block-structured languages. Hence the Pascal version could be useful in producing versions in other structured languages. Rather than use identifiers of minimal length (which therefore have little mnemonic significance), the Pascal version uses meaningful identifiers as follows [Note: A few changes have been made for C]: BASIC C BASIC C BASIC C A J S StickyBit A1 AInverse J0 NoErrors T B Radix [Failure] T0 Underflow B1 BInverse J1 NoErrors T2 ThirtyTwo B2 RadixD2 [SeriousDefect] T5 OneAndHalf B9 BMinusU2 J2 NoErrors T7 TwentySeven C [Defect] T8 TwoForty C1 CInverse J3 NoErrors U OneUlp D [Flaw] U0 UnderflowThreshold D4 FourD K PageNo U1 E0 L Milestone U2 E1 M V E2 Exp2 N V0 E3 N1 V8 E5 MinSqEr O Zero V9 E6 SqEr O1 One W E7 MaxSqEr O2 Two X E8 O3 Three X1 E9 O4 Four X8 F1 MinusOne O5 Five X9 Random1 F2 Half O8 Eight Y F3 Third O9 Nine Y1 F6 P Precision Y2 F9 Q Y9 Random2 G1 GMult Q8 Z G2 GDiv Q9 Z0 PseudoZero G3 GAddSub R Z1 H R1 RMult Z2 H1 HInverse R2 RDiv Z9 I R3 RAddSub IO NoTrials R4 RSqrt I3 IEEE R9 Random9 SqRWrng All the variables in BASIC are true variables and in consequence, the program is more difficult to follow since the "constants" must be determined (the glossary is very helpful). The Pascal version uses Real constants, but checks are added to ensure that the values are correctly converted by the compiler. The major textual change to the Pascal version apart from the identifiersis that named procedures are used, inserting parameters wherehelpful. New procedures are also introduced. The correspondence is as follows: BASIC Pascal lines 90- 140 Pause 170- 250 Instructions 380- 460 Heading 480- 670 Characteristics 690- 870 History 2940-2950 Random 3710-3740 NewD 4040-4080 DoesYequalX 4090-4110 PrintIfNPositive 4640-4850 TestPartialUnderflow */ /* This version of paranoia has been modified to work with GCC's internal software floating point emulation library, as a sanity check of same. I'm doing this in C++ so that I can do operator overloading and not have to modify so damned much of the existing code. */ extern "C" { #include <stdio.h> #include <stddef.h> #include <limits.h> #include <string.h> #include <stdlib.h> #include <math.h> #include <unistd.h> #include <float.h> /* This part is made all the more awful because many gcc headers are not prepared at all to be parsed as C++. The biggest stickler here is const structure members. So we include exactly the pieces that we need. */ #define GTY(x) #include "ansidecl.h" #include "auto-host.h" #include "hwint.h" #undef EXTRA_MODES_FILE struct rtx_def; typedef struct rtx_def *rtx; struct rtvec_def; typedef struct rtvec_def *rtvec; union tree_node; typedef union tree_node *tree; #define DEFTREECODE(SYM, STRING, TYPE, NARGS) SYM, enum tree_code { #include "tree.def" LAST_AND_UNUSED_TREE_CODE }; #undef DEFTREECODE #define ENUM_BITFIELD(X) enum X #define class klass #include "real.h" #undef class } /* We never produce signals from the library. Thus setjmp need do nothing. */ #undef setjmp #define setjmp(x) (0) static bool verbose = false; static int verbose_index = 0; /* ====================================================================== */ /* The implementation of the abstract floating point class based on gcc's real.c. I.e. the object of this exercise. Templated so that we can all fp sizes. */ class real_c_float { public: static const enum machine_mode MODE = SFmode; private: static const int external_max = 128 / 32; static const int internal_max = (sizeof (REAL_VALUE_TYPE) + sizeof (long) + 1) / sizeof (long); long image[external_max < internal_max ? internal_max : external_max]; void from_long(long); void from_str(const char *); void binop(int code, const real_c_float&); void unop(int code); bool cmp(int code, const real_c_float&) const; public: real_c_float() { } real_c_float(long l) { from_long(l); } real_c_float(const char *s) { from_str(s); } real_c_float(const real_c_float &b) { memcpy(image, b.image, sizeof(image)); } const real_c_float& operator= (long l) { from_long(l); return *this; } const real_c_float& operator= (const char *s) { from_str(s); return *this; } const real_c_float& operator= (const real_c_float &b) { memcpy(image, b.image, sizeof(image)); return *this; } const real_c_float& operator+= (const real_c_float &b) { binop(PLUS_EXPR, b); return *this; } const real_c_float& operator-= (const real_c_float &b) { binop(MINUS_EXPR, b); return *this; } const real_c_float& operator*= (const real_c_float &b) { binop(MULT_EXPR, b); return *this; } const real_c_float& operator/= (const real_c_float &b) { binop(RDIV_EXPR, b); return *this; } real_c_float operator- () const { real_c_float r(*this); r.unop(NEGATE_EXPR); return r; } real_c_float abs () const { real_c_float r(*this); r.unop(ABS_EXPR); return r; } bool operator < (const real_c_float &b) const { return cmp(LT_EXPR, b); } bool operator <= (const real_c_float &b) const { return cmp(LE_EXPR, b); } bool operator == (const real_c_float &b) const { return cmp(EQ_EXPR, b); } bool operator != (const real_c_float &b) const { return cmp(NE_EXPR, b); } bool operator >= (const real_c_float &b) const { return cmp(GE_EXPR, b); } bool operator > (const real_c_float &b) const { return cmp(GT_EXPR, b); } const char * str () const; const char * hex () const; long integer () const; int exp () const; void ldexp (int); }; void real_c_float::from_long (long l) { REAL_VALUE_TYPE f; real_from_integer (&f, MODE, l, l < 0 ? -1 : 0, 0); real_to_target (image, &f, MODE); } void real_c_float::from_str (const char *s) { REAL_VALUE_TYPE f; const char *p = s; if (*p == '-' || *p == '+') p++; if (strcasecmp(p, "inf") == 0) { real_inf (&f); if (*s == '-') real_arithmetic (&f, NEGATE_EXPR, &f, NULL); } else if (strcasecmp(p, "nan") == 0) real_nan (&f, "", 1, MODE); else real_from_string (&f, s); real_to_target (image, &f, MODE); } void real_c_float::binop (int code, const real_c_float &b) { REAL_VALUE_TYPE ai, bi, ri; real_from_target (&ai, image, MODE); real_from_target (&bi, b.image, MODE); real_arithmetic (&ri, code, &ai, &bi); real_to_target (image, &ri, MODE); if (verbose) { char ab[64], bb[64], rb[64]; const real_format *fmt = real_format_for_mode[MODE - QFmode]; const int digits = (fmt->p * fmt->log2_b + 3) / 4; char symbol_for_code; real_from_target (&ri, image, MODE); real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0); real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0); real_to_hexadecimal (rb, &ri, sizeof(rb), digits, 0); switch (code) { case PLUS_EXPR: symbol_for_code = '+'; break; case MINUS_EXPR: symbol_for_code = '-'; break; case MULT_EXPR: symbol_for_code = '*'; break; case RDIV_EXPR: symbol_for_code = '/'; break; default: abort (); } fprintf (stderr, "%6d: %s %c %s = %s\n", verbose_index++, ab, symbol_for_code, bb, rb); } } void real_c_float::unop (int code) { REAL_VALUE_TYPE ai, ri; real_from_target (&ai, image, MODE); real_arithmetic (&ri, code, &ai, NULL); real_to_target (image, &ri, MODE); if (verbose) { char ab[64], rb[64]; const real_format *fmt = real_format_for_mode[MODE - QFmode]; const int digits = (fmt->p * fmt->log2_b + 3) / 4; const char *symbol_for_code; real_from_target (&ri, image, MODE); real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0); real_to_hexadecimal (rb, &ri, sizeof(rb), digits, 0); switch (code) { case NEGATE_EXPR: symbol_for_code = "-"; break; case ABS_EXPR: symbol_for_code = "abs "; break; default: abort (); } fprintf (stderr, "%6d: %s%s = %s\n", verbose_index++, symbol_for_code, ab, rb); } } bool real_c_float::cmp (int code, const real_c_float &b) const { REAL_VALUE_TYPE ai, bi; bool ret; real_from_target (&ai, image, MODE); real_from_target (&bi, b.image, MODE); ret = real_compare (code, &ai, &bi); if (verbose) { char ab[64], bb[64]; const real_format *fmt = real_format_for_mode[MODE - QFmode]; const int digits = (fmt->p * fmt->log2_b + 3) / 4; const char *symbol_for_code; real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0); real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0); switch (code) { case LT_EXPR: symbol_for_code = "<"; break; case LE_EXPR: symbol_for_code = "<="; break; case EQ_EXPR: symbol_for_code = "=="; break; case NE_EXPR: symbol_for_code = "!="; break; case GE_EXPR: symbol_for_code = ">="; break; case GT_EXPR: symbol_for_code = ">"; break; default: abort (); } fprintf (stderr, "%6d: %s %s %s = %s\n", verbose_index++, ab, symbol_for_code, bb, (ret ? "true" : "false")); } return ret; } const char * real_c_float::str() const { REAL_VALUE_TYPE f; const real_format *fmt = real_format_for_mode[MODE - QFmode]; const int digits = int(fmt->p * fmt->log2_b * .30102999566398119521 + 1); real_from_target (&f, image, MODE); char *buf = new char[digits + 10]; real_to_decimal (buf, &f, digits+10, digits, 0); return buf; } const char * real_c_float::hex() const { REAL_VALUE_TYPE f; const real_format *fmt = real_format_for_mode[MODE - QFmode]; const int digits = (fmt->p * fmt->log2_b + 3) / 4; real_from_target (&f, image, MODE); char *buf = new char[digits + 10]; real_to_hexadecimal (buf, &f, digits+10, digits, 0); return buf; } long real_c_float::integer() const { REAL_VALUE_TYPE f; real_from_target (&f, image, MODE); return real_to_integer (&f); } int real_c_float::exp() const { REAL_VALUE_TYPE f; real_from_target (&f, image, MODE); return real_exponent (&f); } void real_c_float::ldexp (int exp) { REAL_VALUE_TYPE ai; real_from_target (&ai, image, MODE); real_ldexp (&ai, &ai, exp); real_to_target (image, &ai, MODE); } /* ====================================================================== */ /* An implementation of the abstract floating point class that uses native arithmetic. Exists for reference and debugging. */ template<typename T> class native_float { private: // Force intermediate results back to memory. volatile T image; static T from_str (const char *); static T do_abs (T); static T verbose_binop (T, char, T, T); static T verbose_unop (const char *, T, T); static bool verbose_cmp (T, const char *, T, bool); public: native_float() { } native_float(long l) { image = l; } native_float(const char *s) { image = from_str(s); } native_float(const native_float &b) { image = b.image; } const native_float& operator= (long l) { image = l; return *this; } const native_float& operator= (const char *s) { image = from_str(s); return *this; } const native_float& operator= (const native_float &b) { image = b.image; return *this; } const native_float& operator+= (const native_float &b) { image = verbose_binop(image, '+', b.image, image + b.image); return *this; } const native_float& operator-= (const native_float &b) { image = verbose_binop(image, '-', b.image, image - b.image); return *this; } const native_float& operator*= (const native_float &b) { image = verbose_binop(image, '*', b.image, image * b.image); return *this; } const native_float& operator/= (const native_float &b) { image = verbose_binop(image, '/', b.image, image / b.image); return *this; } native_float operator- () const { native_float r; r.image = verbose_unop("-", image, -image); return r; } native_float abs () const { native_float r; r.image = verbose_unop("abs ", image, do_abs(image)); return r; } bool operator < (const native_float &b) const { return verbose_cmp(image, "<", b.image, image < b.image); } bool operator <= (const native_float &b) const { return verbose_cmp(image, "<=", b.image, image <= b.image); } bool operator == (const native_float &b) const { return verbose_cmp(image, "==", b.image, image == b.image); } bool operator != (const native_float &b) const { return verbose_cmp(image, "!=", b.image, image != b.image); } bool operator >= (const native_float &b) const { return verbose_cmp(image, ">=", b.image, image >= b.image); } bool operator > (const native_float &b) const { return verbose_cmp(image, ">", b.image, image > b.image); } const char * str () const; const char * hex () const; long integer () const { return long(image); } int exp () const; void ldexp (int); }; template<typename T> inline T native_float<T>::from_str (const char *s) { return strtold (s, NULL); } template<> inline float native_float<float>::from_str (const char *s) { return strtof (s, NULL); } template<> inline double native_float<double>::from_str (const char *s) { return strtod (s, NULL); } template<typename T> inline T native_float<T>::do_abs (T image) { return fabsl (image); } template<> inline float native_float<float>::do_abs (float image) { return fabsf (image); } template<> inline double native_float<double>::do_abs (double image) { return fabs (image); } template<typename T> T native_float<T>::verbose_binop (T a, char symbol, T b, T r) { if (verbose) { const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1; #ifdef NO_LONG_DOUBLE fprintf (stderr, "%6d: %.*a %c %.*a = %.*a\n", verbose_index++, digits, (double)a, symbol, digits, (double)b, digits, (double)r); #else fprintf (stderr, "%6d: %.*La %c %.*La = %.*La\n", verbose_index++, digits, (long double)a, symbol, digits, (long double)b, digits, (long double)r); #endif } return r; } template<typename T> T native_float<T>::verbose_unop (const char *symbol, T a, T r) { if (verbose) { const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1; #ifdef NO_LONG_DOUBLE fprintf (stderr, "%6d: %s%.*a = %.*a\n", verbose_index++, symbol, digits, (double)a, digits, (double)r); #else fprintf (stderr, "%6d: %s%.*La = %.*La\n", verbose_index++, symbol, digits, (long double)a, digits, (long double)r); #endif } return r; } template<typename T> bool native_float<T>::verbose_cmp (T a, const char *symbol, T b, bool r) { if (verbose) { const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1; #ifndef NO_LONG_DOUBLE fprintf (stderr, "%6d: %.*a %s %.*a = %s\n", verbose_index++, digits, (double)a, symbol, digits, (double)b, (r ? "true" : "false")); #else fprintf (stderr, "%6d: %.*La %s %.*La = %s\n", verbose_index++, digits, (long double)a, symbol, digits, (long double)b, (r ? "true" : "false")); #endif } return r; } template<typename T> const char * native_float<T>::str() const { char *buf = new char[50]; const int digits = int(sizeof(T) * CHAR_BIT * .30102999566398119521 + 1); #ifndef NO_LONG_DOUBLE sprintf (buf, "%.*e", digits - 1, (double) image); #else sprintf (buf, "%.*Le", digits - 1, (long double) image); #endif return buf; } template<typename T> const char * native_float<T>::hex() const { char *buf = new char[50]; const int digits = int(sizeof(T) * CHAR_BIT / 4); #ifndef NO_LONG_DOUBLE sprintf (buf, "%.*a", digits - 1, (double) image); #else sprintf (buf, "%.*La", digits - 1, (long double) image); #endif return buf; } template<typename T> int native_float<T>::exp() const { int e; frexp (image, &e); return e; } template<typename T> void native_float<T>::ldexp (int exp) { image = ldexpl (image, exp); } template<> void native_float<float>::ldexp (int exp) { image = ldexpf (image, exp); } template<> void native_float<double>::ldexp (int exp) { image = ::ldexp (image, exp); } /* ====================================================================== */ /* Some libm routines that Paranoia expects to be available. */ template<typename FLOAT> inline FLOAT FABS (const FLOAT &f) { return f.abs(); } template<typename FLOAT, typename RHS> inline FLOAT operator+ (const FLOAT &a, const RHS &b) { return FLOAT(a) += FLOAT(b); } template<typename FLOAT, typename RHS> inline FLOAT operator- (const FLOAT &a, const RHS &b) { return FLOAT(a) -= FLOAT(b); } template<typename FLOAT, typename RHS> inline FLOAT operator* (const FLOAT &a, const RHS &b) { return FLOAT(a) *= FLOAT(b); } template<typename FLOAT, typename RHS> inline FLOAT operator/ (const FLOAT &a, const RHS &b) { return FLOAT(a) /= FLOAT(b); } template<typename FLOAT> FLOAT FLOOR (const FLOAT &f) { /* ??? This is only correct when F is representable as an integer. */ long i = f.integer(); FLOAT r; r = i; if (i < 0 && f != r) r = i - 1; return r; } template<typename FLOAT> FLOAT SQRT (const FLOAT &f) { #if 0 FLOAT zero = long(0); FLOAT two = 2; FLOAT one = 1; FLOAT diff, diff2; FLOAT z, t; if (f == zero) return zero; if (f < zero) return zero / zero; if (f == one) return f; z = f; z.ldexp (-f.exp() / 2); diff2 = FABS (z * z - f); if (diff2 > zero) while (1) { t = (f / (two * z)) + (z / two); diff = FABS (t * t - f); if (diff >= diff2) break; z = t; diff2 = diff; } return z; #elif defined(NO_LONG_DOUBLE) double d; char buf[64]; d = strtod (f.hex(), NULL); d = sqrt (d); sprintf(buf, "%.35a", d); return FLOAT(buf); #else long double ld; char buf[64]; ld = strtold (f.hex(), NULL); ld = sqrtl (ld); sprintf(buf, "%.35La", ld); return FLOAT(buf); #endif } template<typename FLOAT> FLOAT LOG (FLOAT x) { #if 0 FLOAT zero = long(0); FLOAT one = 1; if (x <= zero) return zero / zero; if (x == one) return zero; int exp = x.exp() - 1; x.ldexp(-exp); FLOAT xm1 = x - one; FLOAT y = xm1; long n = 2; FLOAT sum = xm1; while (1) { y *= xm1; FLOAT term = y / FLOAT (n); FLOAT next = sum + term; if (next == sum) break; sum = next; if (++n == 1000) break; } if (exp) sum += FLOAT (exp) * FLOAT(".69314718055994530941"); return sum; #elif defined (NO_LONG_DOUBLE) double d; char buf[64]; d = strtod (x.hex(), NULL); d = log (d); sprintf(buf, "%.35a", d); return FLOAT(buf); #else long double ld; char buf[64]; ld = strtold (x.hex(), NULL); ld = logl (ld); sprintf(buf, "%.35La", ld); return FLOAT(buf); #endif } template<typename FLOAT> FLOAT EXP (const FLOAT &x) { /* Cheat. */ #ifdef NO_LONG_DOUBLE double d; char buf[64]; d = strtod (x.hex(), NULL); d = exp (d); sprintf(buf, "%.35a", d); return FLOAT(buf); #else long double ld; char buf[64]; ld = strtold (x.hex(), NULL); ld = expl (ld); sprintf(buf, "%.35La", ld); return FLOAT(buf); #endif } template<typename FLOAT> FLOAT POW (const FLOAT &base, const FLOAT &exp) { /* Cheat. */ #ifdef NO_LONG_DOUBLE double d1, d2; char buf[64]; d1 = strtod (base.hex(), NULL); d2 = strtod (exp.hex(), NULL); d1 = pow (d1, d2); sprintf(buf, "%.35a", d1); return FLOAT(buf); #else long double ld1, ld2; char buf[64]; ld1 = strtold (base.hex(), NULL); ld2 = strtold (exp.hex(), NULL); ld1 = powl (ld1, ld2); sprintf(buf, "%.35La", ld1); return FLOAT(buf); #endif } /* ====================================================================== */ /* Real Paranoia begins again here. We wrap the thing in a template so that we can instantiate it for each floating point type we care for. */ int NoTrials = 20; /*Number of tests for commutativity. */ bool do_pause = false; enum Guard { No, Yes }; enum Rounding { Other, Rounded, Chopped }; enum Class { Failure, Serious, Defect, Flaw }; template<typename FLOAT> struct Paranoia { FLOAT Radix, BInvrse, RadixD2, BMinusU2; /* Small floating point constants. */ FLOAT Zero; FLOAT Half; FLOAT One; FLOAT Two; FLOAT Three; FLOAT Four; FLOAT Five; FLOAT Eight; FLOAT Nine; FLOAT TwentySeven; FLOAT ThirtyTwo; FLOAT TwoForty; FLOAT MinusOne; FLOAT OneAndHalf; /* Declarations of Variables. */ int Indx; char ch[8]; FLOAT AInvrse, A1; FLOAT C, CInvrse; FLOAT D, FourD; FLOAT E0, E1, Exp2, E3, MinSqEr; FLOAT SqEr, MaxSqEr, E9; FLOAT Third; FLOAT F6, F9; FLOAT H, HInvrse; int I; FLOAT StickyBit, J; FLOAT MyZero; FLOAT Precision; FLOAT Q, Q9; FLOAT R, Random9; FLOAT T, Underflow, S; FLOAT OneUlp, UfThold, U1, U2; FLOAT V, V0, V9; FLOAT W; FLOAT X, X1, X2, X8, Random1; FLOAT Y, Y1, Y2, Random2; FLOAT Z, PseudoZero, Z1, Z2, Z9; int ErrCnt[4]; int Milestone; int PageNo; int M, N, N1; Guard GMult, GDiv, GAddSub; Rounding RMult, RDiv, RAddSub, RSqrt; int Break, Done, NotMonot, Monot, Anomaly, IEEE, SqRWrng, UfNGrad; /* 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 */ int main (); FLOAT Sign (FLOAT); FLOAT Random (); void Pause (); void BadCond (int, const char *); void SqXMinX (int); void TstCond (int, int, const char *); void notify (const char *); void IsYeqX (); void NewD (); void PrintIfNPositive (); void SR3750 (); void TstPtUf (); // Pretend we're bss. Paranoia() { memset(this, 0, sizeof (*this)); } }; template<typename FLOAT> int Paranoia<FLOAT>::main() { /* First two assignments use integer right-hand sides. */ Zero = long(0); One = long(1); Two = long(2); Three = long(3); Four = long(4); Five = long(5); Eight = long(8); Nine = long(9); TwentySeven = long(27); ThirtyTwo = long(32); TwoForty = long(240); MinusOne = long(-1); Half = "0x1p-1"; OneAndHalf = "0x3p-1"; ErrCnt[Failure] = 0; ErrCnt[Serious] = 0; ErrCnt[Defect] = 0; ErrCnt[Flaw] = 0; PageNo = 1; /*=============================================*/ Milestone = 7; /*=============================================*/ printf ("Program is now RUNNING tests on small integers:\n"); TstCond (Failure, (Zero + Zero == Zero), "0+0 != 0"); TstCond (Failure, (One - One == Zero), "1-1 != 0"); TstCond (Failure, (One > Zero), "1 <= 0"); TstCond (Failure, (One + One == Two), "1+1 != 2"); Z = -Zero; if (Z != Zero) { ErrCnt[Failure] = ErrCnt[Failure] + 1; printf ("Comparison alleges that -0.0 is Non-zero!\n"); U2 = "0.001"; Radix = 1; TstPtUf (); } TstCond (Failure, (Three == Two + One), "3 != 2+1"); TstCond (Failure, (Four == Three + One), "4 != 3+1"); TstCond (Failure, (Four + Two * (-Two) == Zero), "4 + 2*(-2) != 0"); TstCond (Failure, (Four - Three - One == Zero), "4-3-1 != 0"); TstCond (Failure, (MinusOne == (Zero - One)), "-1 != 0-1"); TstCond (Failure, (MinusOne + One == Zero), "-1+1 != 0"); TstCond (Failure, (One + MinusOne == Zero), "1+(-1) != 0"); TstCond (Failure, (MinusOne + FABS (One) == Zero), "-1+abs(1) != 0"); TstCond (Failure, (MinusOne + MinusOne * MinusOne == Zero), "-1+(-1)*(-1) != 0"); TstCond (Failure, Half + MinusOne + Half == Zero, "1/2 + (-1) + 1/2 != 0"); /*=============================================*/ Milestone = 10; /*=============================================*/ TstCond (Failure, (Nine == Three * Three), "9 != 3*3"); TstCond (Failure, (TwentySeven == Nine * Three), "27 != 9*3"); TstCond (Failure, (Eight == Four + Four), "8 != 4+4"); TstCond (Failure, (ThirtyTwo == Eight * Four), "32 != 8*4"); TstCond (Failure, (ThirtyTwo - TwentySeven - Four - One == Zero), "32-27-4-1 != 0"); TstCond (Failure, Five == Four + One, "5 != 4+1"); TstCond (Failure, TwoForty == Four * Five * Three * Four, "240 != 4*5*3*4"); TstCond (Failure, TwoForty / Three - Four * Four * Five == Zero, "240/3 - 4*4*5 != 0"); TstCond (Failure, TwoForty / Four - Five * Three * Four == Zero, "240/4 - 5*3*4 != 0"); TstCond (Failure, TwoForty / Five - Four * Three * Four == Zero, "240/5 - 4*3*4 != 0"); if (ErrCnt[Failure] == 0) { printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n"); printf ("\n"); } printf ("Searching for Radix and Precision.\n"); W = One; do { W = W + W; Y = W + One; Z = Y - W; Y = Z - One; } while (MinusOne + FABS (Y) < Zero); /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ... */ Precision = Zero; Y = One; do { Radix = W + Y; Y = Y + Y; Radix = Radix - W; } while (Radix == Zero); if (Radix < Two) Radix = One; printf ("Radix = %s .\n", Radix.str()); if (Radix != One) { W = One; do { Precision = Precision + One; W = W * Radix; Y = W + One; } while ((Y - W) == One); } /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 ... */ U1 = One / W; U2 = Radix * U1; printf ("Closest relative separation found is U1 = %s .\n\n", U1.str()); printf ("Recalculating radix and precision\n "); /*save old values */ E0 = Radix; E1 = U1; E9 = U2; E3 = Precision; X = Four / Three; Third = X - One; F6 = Half - Third; X = F6 + F6; X = FABS (X - Third); if (X < U2) X = U2; /*... now X = (unknown no.) ulps of 1+... */ do { U2 = X; Y = Half * U2 + ThirtyTwo * U2 * U2; Y = One + Y; X = Y - One; } while (!((U2 <= X) || (X <= Zero))); /*... now U2 == 1 ulp of 1 + ... */ X = Two / Three; F6 = X - Half; Third = F6 + F6; X = Third - Half; X = FABS (X + F6); if (X < U1) X = U1; /*... now X == (unknown no.) ulps of 1 -... */ do { U1 = X; Y = Half * U1 + ThirtyTwo * U1 * U1; Y = Half - Y; X = Half + Y; Y = Half - X; X = Half + Y; } while (!((U1 <= X) || (X <= Zero))); /*... now U1 == 1 ulp of 1 - ... */ if (U1 == E1) printf ("confirms closest relative separation U1 .\n"); else printf ("gets better closest relative separation U1 = %s .\n", U1.str()); W = One / U1; F9 = (Half - U1) + Half; Radix = FLOOR (FLOAT ("0.01") + U2 / U1); if (Radix == E0) printf ("Radix confirmed.\n"); else printf ("MYSTERY: recalculated Radix = %s .\n", Radix.str()); TstCond (Defect, Radix <= Eight + Eight, "Radix is too big: roundoff problems"); TstCond (Flaw, (Radix == Two) || (Radix == 10) || (Radix == One), "Radix is not as good as 2 or 10"); /*=============================================*/ Milestone = 20; /*=============================================*/ TstCond (Failure, F9 - Half < Half, "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?"); X = F9; I = 1; Y = X - Half; Z = Y - Half; TstCond (Failure, (X != One) || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0"); X = One + U2; I = 0; /*=============================================*/ Milestone = 25; /*=============================================*/ /*... BMinusU2 = nextafter(Radix, 0) */ BMinusU2 = Radix - One; BMinusU2 = (BMinusU2 - U2) + One; /* Purify Integers */ if (Radix != One) { X = -TwoForty * LOG (U1) / LOG (Radix); Y = FLOOR (Half + X); if (FABS (X - Y) * Four < One) X = Y; Precision = X / TwoForty; Y = FLOOR (Half + Precision); if (FABS (Precision - Y) * TwoForty < Half) Precision = Y; } if ((Precision != FLOOR (Precision)) || (Radix == One)) { printf ("Precision cannot be characterized by an Integer number\n"); printf ("of significant digits but, by itself, this is a minor flaw.\n"); } if (Radix == One) printf ("logarithmic encoding has precision characterized solely by U1.\n"); else printf ("The number of significant digits of the Radix is %s .\n", Precision.str()); TstCond (Serious, U2 * Nine * Nine * TwoForty < One, "Precision worse than 5 decimal figures "); /*=============================================*/ Milestone = 30; /*=============================================*/ /* Test for extra-precise subexpressions */ X = FABS (((Four / Three - One) - One / Four) * Three - One / Four); do { Z2 = X; X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One; } while (!((Z2 <= X) || (X <= Zero))); X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four); do { Z1 = Z; Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1)) + One / Two)) + One / Two; } while (!((Z1 <= Z) || (Z <= Zero))); do { do { Y1 = Y; Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half)) + Half; } while (!((Y1 <= Y) || (Y <= Zero))); X1 = X; X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9; } while (!((X1 <= X) || (X <= Zero))); if ((X1 != Y1) || (X1 != Z1)) { BadCond (Serious, "Disagreements among the values X1, Y1, Z1,\n"); printf ("respectively %s, %s, %s,\n", X1.str(), Y1.str(), Z1.str()); printf ("are symptoms of inconsistencies introduced\n"); printf ("by extra-precise evaluation of arithmetic subexpressions.\n"); notify ("Possibly some part of this"); if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf ("That feature is not tested further by this program.\n"); } else { if ((Z1 != U1) || (Z2 != U2)) { if ((Z1 >= U1) || (Z2 >= U2)) { BadCond (Failure, ""); notify ("Precision"); printf ("\tU1 = %s, Z1 - U1 = %s\n", U1.str(), (Z1 - U1).str()); printf ("\tU2 = %s, Z2 - U2 = %s\n", U2.str(), (Z2 - U2).str()); } else { if ((Z1 <= Zero) || (Z2 <= Zero)) { printf ("Because of unusual Radix = %s", Radix.str()); printf (", or exact rational arithmetic a result\n"); printf ("Z1 = %s, or Z2 = %s ", Z1.str(), Z2.str()); notify ("of an\nextra-precision"); } if (Z1 != Z2 || Z1 > Zero) { X = Z1 / U1; Y = Z2 / U2; if (Y > X) X = Y; Q = -LOG (X); printf ("Some subexpressions appear to be calculated " "extra precisely\n"); printf ("with about %s extra B-digits, i.e.\n", (Q / LOG (Radix)).str()); printf ("roughly %s extra significant decimals.\n", (Q / LOG (FLOAT (10))).str()); } printf ("That feature is not tested further by this program.\n"); } } } Pause (); /*=============================================*/ Milestone = 35; /*=============================================*/ if (Radix >= Two) { X = W / (Radix * Radix); Y = X + One; Z = Y - X; T = Z + U2; X = T - Z; TstCond (Failure, X == U2, "Subtraction is not normalized X=Y,X+Z != Y+Z!"); if (X == U2) printf ("Subtraction appears to be normalized, as it should be."); } printf ("\nChecking for guard digit in *, /, and -.\n"); Y = F9 * One; Z = One * F9; X = F9 - Half; Y = (Y - Half) - X; Z = (Z - Half) - X; X = One + U2; T = X * Radix; R = Radix * X; X = T - Radix; X = X - Radix * U2; T = R - Radix; T = T - Radix * U2; X = X * (Radix - One); T = T * (Radix - One); if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes; else { GMult = No; TstCond (Serious, false, "* lacks a Guard Digit, so 1*X != X"); } Z = Radix * U2; X = One + Z; Y = FABS ((X + Z) - X * X) - U2; X = One - U2; Z = FABS ((X - U2) - X * X) - U1; TstCond (Failure, (Y <= Zero) && (Z <= Zero), "* gets too many final digits wrong.\n"); Y = One - U2; X = One + U2; Z = One / Y; Y = Z - X; X = One / Three; Z = Three / Nine; X = X - Z; T = Nine / TwentySeven; Z = Z - T; TstCond (Defect, X == Zero && Y == Zero && Z == Zero, "Division lacks a Guard Digit, so error can exceed 1 ulp\n" "or 1/3 and 3/9 and 9/27 may disagree"); Y = F9 / One; X = F9 - Half; Y = (Y - Half) - X; X = One + U2; T = X / One; X = T - X; if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes; else { GDiv = No; TstCond (Serious, false, "Division lacks a Guard Digit, so X/1 != X"); } X = One / (One + U2); Y = X - Half - Half; TstCond (Serious, Y < Zero, "Computed value of 1/1.000..1 >= 1"); X = One - U2; Y = One + Radix * U2; Z = X * Radix; T = Y * Radix; R = Z / Radix; StickyBit = T / Radix; X = R - X; Y = StickyBit - Y; TstCond (Failure, X == Zero && Y == Zero, "* and/or / gets too many last digits wrong"); Y = One - U1; X = One - F9; Y = One - Y; T = Radix - U2; Z = Radix - BMinusU2; T = Radix - T; if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes; else { GAddSub = No; TstCond (Serious, false, "- lacks Guard Digit, so cancellation is obscured"); } if (F9 != One && F9 - One >= Zero) { BadCond (Serious, "comparison alleges (1-U1) < 1 although\n"); printf (" subtraction 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; /*=============================================*/ Pause (); printf ("Checking rounding on multiply, divide and add/subtract.\n"); RMult = Other; RDiv = Other; RAddSub = Other; RadixD2 = Radix / Two; A1 = Two; Done = false; do { AInvrse = Radix; do { X = AInvrse; AInvrse = AInvrse / A1; } while (!(FLOOR (AInvrse) != AInvrse)); Done = (X == One) || (A1 > Three); if (!Done) A1 = Nine + One; } while (!(Done)); if (X == One) A1 = Radix; AInvrse = One / A1; X = A1; Y = AInvrse; Done = false; do { Z = X * Y - Half; TstCond (Failure, Z == Half, "X * (1/X) differs from 1"); Done = X == Radix; X = Radix; Y = One / X; } while (!(Done)); Y2 = One + U2; Y1 = One - U2; X = OneAndHalf - U2; Y = OneAndHalf + U2; Z = (X - U2) * Y2; T = Y * Y1; Z = Z - X; T = T - X; X = X * Y2; Y = (Y + U2) * Y1; X = X - OneAndHalf; Y = Y - OneAndHalf; if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) { X = (OneAndHalf + U2) * Y2; Y = OneAndHalf - U2 - U2; Z = OneAndHalf + U2 + U2; T = (OneAndHalf - U2) * Y1; X = X - (Z + U2); StickyBit = Y * Y1; S = Z * Y2; T = T - Y; Y = (U2 - Y) + StickyBit; Z = S - (Z + U2 + U2); StickyBit = (Y2 + U2) * Y1; Y1 = Y2 * Y1; StickyBit = StickyBit - Y2; Y1 = Y1 - Half; if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero) && (StickyBit == Zero) && (Y1 == Half)) { RMult = Rounded; printf ("Multiplication appears to round correctly.\n"); } else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero) && (T < Zero) && (StickyBit + U2 == Zero) && (Y1 < Half)) { RMult = Chopped; printf ("Multiplication appears to chop.\n"); } else printf ("* is neither chopped nor correctly rounded.\n"); if ((RMult == Rounded) && (GMult == No)) notify ("Multiplication"); } else printf ("* is neither chopped nor correctly rounded.\n"); /*=============================================*/ Milestone = 45; /*=============================================*/ Y2 = One + U2; Y1 = One - U2; Z = OneAndHalf + U2 + U2; X = Z / Y2; T = OneAndHalf - U2 - U2; Y = (T - U2) / Y1; Z = (Z + U2) / Y2; X = X - OneAndHalf; Y = Y - T; T = T / Y1; Z = Z - (OneAndHalf + U2); T = (U2 - OneAndHalf) + T; if (!((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) { X = OneAndHalf / Y2; Y = OneAndHalf - U2; Z = OneAndHalf + U2; X = X - Y; T = OneAndHalf / Y1; Y = Y / Y1; T = T - (Z + U2); Y = Y - Z; Z = Z / Y2; Y1 = (Y2 + U2) / Y2; Z = Z - OneAndHalf; Y2 = Y1 - Y2; Y1 = (F9 - U1) / F9; if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero) && (Y2 == Zero) && (Y2 == Zero) && (Y1 - Half == F9 - Half)) { RDiv = Rounded; printf ("Division appears to round correctly.\n"); if (GDiv == No) notify ("Division"); } else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero) && (Y2 < Zero) && (Y1 - Half < F9 - Half)) { RDiv = Chopped; printf ("Division appears to chop.\n"); } } if (RDiv == Other) printf ("/ is neither chopped nor correctly rounded.\n"); BInvrse = One / Radix; TstCond (Failure, (BInvrse * Radix - Half == Half), "Radix * ( 1 / Radix ) differs from 1"); /*=============================================*/ Milestone = 50; /*=============================================*/ TstCond (Failure, ((F9 + U1) - Half == Half) && ((BMinusU2 + U2) - One == Radix - One), "Incomplete carry-propagation in Addition"); X = One - U1 * U1; Y = One + U2 * (One - U2); Z = F9 - Half; X = (X - Half) - Z; Y = Y - One; if ((X == Zero) && (Y == Zero)) { RAddSub = Chopped; printf ("Add/Subtract appears to be chopped.\n"); } if (GAddSub == Yes) { X = (Half + U2) * U2; Y = (Half - U2) * U2; X = One + X; Y = One + Y; X = (One + U2) - X; Y = One - Y; if ((X == Zero) && (Y == Zero)) { X = (Half + U2) * U1; Y = (Half - U2) * U1; X = One - X; Y = One - Y; X = F9 - X; Y = One - Y; if ((X == Zero) && (Y == Zero)) { RAddSub = Rounded; printf ("Addition/Subtraction appears to round correctly.\n"); if (GAddSub == No) notify ("Add/Subtract"); } 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"); S = One; X = One + Half * (One + Half); Y = (One + U2) * Half; Z = X - Y; T = Y - X; StickyBit = Z + T; if (StickyBit != Zero) { S = Zero; BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n"); } StickyBit = Zero; if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes) && (RMult == Rounded) && (RDiv == Rounded) && (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2)) { printf ("Checking for sticky bit.\n"); X = (Half + U1) * U2; Y = Half * U2; Z = One + Y; T = One + X; if ((Z - One <= Zero) && (T - One >= U2)) { Z = T + Y; Y = Z - X; if ((Z - T >= U2) && (Y - T == Zero)) { X = (Half + U1) * U1; Y = Half * U1; Z = One - Y; T = One - X; if ((Z - One == Zero) && (T - F9 == Zero)) { Z = (Half - U1) * U1; T = F9 - Z; Q = F9 - Y; if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) { Z = (One + U2) * OneAndHalf; T = (OneAndHalf + U2) - Z + U2; X = One + Half / Radix; Y = One + Radix * U2; Z = X * Y; if (T == Zero && X + Radix * U2 - Z == Zero) { if (Radix != Two) { X = Two + U2; Y = X / Two; if ((Y - One == Zero)) StickyBit = S; } else StickyBit = S; } } } } } } if (StickyBit == One) printf ("Sticky bit apparently used correctly.\n"); else printf ("Sticky bit used incorrectly or not at all.\n"); TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No || RMult == Other || RDiv == Other || RAddSub == Other), "lack(s) of guard digits or failure(s) to correctly round or chop\n\ (noted above) count as one flaw in the final tally below"); /*=============================================*/ Milestone = 60; /*=============================================*/ printf ("\n"); printf ("Does Multiplication commute? "); printf ("Testing on %d random pairs.\n", NoTrials); Random9 = SQRT (FLOAT (3)); Random1 = Third; I = 1; do { X = Random (); Y = Random (); Z9 = Y * X; Z = X * Y; Z9 = Z - Z9; I = I + 1; } while (!((I > NoTrials) || (Z9 != Zero))); if (I == NoTrials) { Random1 = One + Half / Three; Random2 = (U2 + U1) + One; Z = Random1 * Random2; Y = Random2 * Random1; Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half / Three) * ((U2 + U1) + One); } if (!((I == NoTrials) || (Z9 == Zero))) BadCond (Defect, "X * Y == Y * X trial fails.\n"); else printf (" No failures found in %d integer pairs.\n", NoTrials); /*=============================================*/ Milestone = 70; /*=============================================*/ printf ("\nRunning test of square root(x).\n"); TstCond (Failure, (Zero == SQRT (Zero)) && (-Zero == SQRT (-Zero)) && (One == SQRT (One)), "Square root of 0.0, -0.0 or 1.0 wrong"); MinSqEr = Zero; MaxSqEr = Zero; J = Zero; X = Radix; OneUlp = U2; SqXMinX (Serious); X = BInvrse; OneUlp = BInvrse * U1; SqXMinX (Serious); X = U1; OneUlp = U1 * U1; SqXMinX (Serious); if (J != Zero) Pause (); printf ("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials); J = Zero; X = Two; Y = Radix; if ((Radix != One)) do { X = Y; Y = Radix * Y; } while (!((Y - X >= NoTrials))); OneUlp = X * U2; I = 1; while (I <= NoTrials) { X = X + One; SqXMinX (Defect); if (J > Zero) break; I = I + 1; } printf ("Test for sqrt monotonicity.\n"); I = -1; X = BMinusU2; Y = Radix; Z = Radix + Radix * U2; NotMonot = false; Monot = false; while (!(NotMonot || Monot)) { I = I + 1; X = SQRT (X); Q = SQRT (Y); Z = SQRT (Z); if ((X > Q) || (Q > Z)) NotMonot = true; else { Q = FLOOR (Q + Half); if (!(I > 0 || Radix == Q * Q)) Monot = true; else if (I > 0) { if (I > 1) Monot = true; else { Y = Y * BInvrse; X = Y - U1; Z = Y + U1; } } else { Y = Q; X = Y - U2; Z = Y + U2; } } } if (Monot) printf ("sqrt has passed a test for Monotonicity.\n"); else { BadCond (Defect, ""); printf ("sqrt(X) is non-monotonic for X near %s .\n", Y.str()); } /*=============================================*/ Milestone = 110; /*=============================================*/ printf ("Seeking Underflow thresholds UfThold and E0.\n"); D = U1; if (Precision != FLOOR (Precision)) { D = BInvrse; X = Precision; do { D = D * BInvrse; X = X - One; } while (X > Zero); } Y = One; Z = D; /* ... D is power of 1/Radix < 1. */ do { C = Y; Y = Z; Z = Y * Y; } while ((Y > Z) && (Z + Z > Z)); Y = C; Z = Y * D; do { C = Y; Y = Z; Z = Y * D; } while ((Y > Z) && (Z + Z > Z)); if (Radix < Two) HInvrse = Two; else HInvrse = Radix; H = One / HInvrse; /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */ CInvrse = One / C; E0 = C; Z = E0 * H; /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */ do { Y = E0; E0 = Z; Z = E0 * H; } while ((E0 > Z) && (Z + Z > Z)); UfThold = E0; E1 = Zero; Q = Zero; E9 = U2; S = One + E9; D = C * S; if (D <= C) { E9 = Radix * U2; S = One + E9; D = C * S; if (D <= C) { BadCond (Failure, "multiplication gets too many last digits wrong.\n"); Underflow = E0; Y1 = Zero; PseudoZero = Z; Pause (); } } else { Underflow = D; PseudoZero = Underflow * H; UfThold = Zero; do { Y1 = Underflow; Underflow = PseudoZero; if (E1 + E1 <= E1) { Y2 = Underflow * HInvrse; E1 = FABS (Y1 - Y2); Q = Y1; if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1; } PseudoZero = PseudoZero * H; } while ((Underflow > PseudoZero) && (PseudoZero + PseudoZero > PseudoZero)); } /* Comment line 4530 .. 4560 */ if (PseudoZero != Zero) { printf ("\n"); Z = PseudoZero; /* ... Test PseudoZero for "phoney- zero" violates */ /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero ... */ if (PseudoZero <= Zero) { BadCond (Failure, "Positive expressions can underflow to an\n"); printf ("allegedly negative value\n"); printf ("PseudoZero that prints out as: %s .\n", PseudoZero.str()); X = -PseudoZero; if (X <= Zero) { printf ("But -PseudoZero, which should be\n"); printf ("positive, isn't; it prints out as %s .\n", X.str()); } } else { BadCond (Flaw, "Underflow can stick at an allegedly positive\n"); printf ("value PseudoZero that prints out as %s .\n", PseudoZero.str()); } TstPtUf (); } /*=============================================*/ Milestone = 120; /*=============================================*/ if (CInvrse * Y > CInvrse * Y1) { S = H * S; E0 = Underflow; } if (!((E1 == Zero) || (E1 == E0))) { BadCond (Defect, ""); if (E1 < E0) { printf ("Products underflow at a higher"); printf (" threshold than differences.\n"); if (PseudoZero == Zero) E0 = E1; } else { printf ("Difference underflows at a higher"); printf (" threshold than products.\n"); } } printf ("Smallest strictly positive number found is E0 = %s .\n", E0.str()); Z = E0; TstPtUf (); Underflow = E0; if (N == 1) Underflow = Y; I = 4; if (E1 == Zero) I = 3; if (UfThold == Zero) I = I - 2; UfNGrad = true; switch (I) { case 1: UfThold = Underflow; if ((CInvrse * Q) != ((CInvrse * Y) * S)) { UfThold = Y; BadCond (Failure, "Either accuracy deteriorates as numbers\n"); printf ("approach a threshold = %s\n", UfThold.str()); printf (" coming down from %s\n", C.str()); printf (" or else multiplication gets too many last digits wrong.\n"); } Pause (); break; case 2: BadCond (Failure, "Underflow confuses Comparison, which alleges that\n"); printf ("Q == Y while denying that |Q - Y| == 0; these values\n"); printf ("print out as Q = %s, Y = %s .\n", Q.str(), Y2.str()); printf ("|Q - Y| = %s .\n", FABS (Q - Y2).str()); UfThold = Q; break; case 3: X = X; break; case 4: if ((Q == UfThold) && (E1 == E0) && (FABS (UfThold - E1 / E9) <= E1)) { UfNGrad = false; printf ("Underflow is gradual; it incurs Absolute Error =\n"); printf ("(roundoff in UfThold) < E0.\n"); Y = E0 * CInvrse; Y = Y * (OneAndHalf + U2); X = CInvrse * (One + U2); Y = Y / X; IEEE = (Y == E0); } } if (UfNGrad) { printf ("\n"); if (setjmp (ovfl_buf)) { printf ("Underflow / UfThold failed!\n"); R = H + H; } else R = SQRT (Underflow / UfThold); if (R <= H) { Z = R * UfThold; X = Z * (One + R * H * (One + H)); } else { Z = UfThold; X = Z * (One + H * H * (One + H)); } if (!((X == Z) || (X - Z != Zero))) { BadCond (Flaw, ""); printf ("X = %s\n\tis not equal to Z = %s .\n", X.str(), Z.str()); Z9 = X - Z; printf ("yet X - Z yields %s .\n", Z9.str()); 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"); if (setjmp (ovfl_buf)) printf ("X / Z fails!\n"); else printf ("X / Z = 1 + %s .\n", ((X / Z - Half) - Half).str()); } } printf ("The Underflow threshold is %s, below which\n", UfThold.str()); printf ("calculation may suffer larger Relative error than "); printf ("merely roundoff.\n"); Y2 = U1 * U1; Y = Y2 * Y2; Y2 = Y * U1; if (Y2 <= UfThold) { if (Y > E0) { BadCond (Defect, ""); I = 5; } else { BadCond (Serious, ""); I = 4; } printf ("Range is too narrow; U1^%d Underflows.\n", I); } /*=============================================*/ Milestone = 130; /*=============================================*/ Y = -FLOOR (Half - TwoForty * LOG (UfThold) / LOG (HInvrse)) / TwoForty; Y2 = Y + Y; printf ("Since underflow occurs below the threshold\n"); printf ("UfThold = (%s) ^ (%s)\nonly underflow ", HInvrse.str(), Y.str()); printf ("should afflict the expression\n\t(%s) ^ (%s);\n", HInvrse.str(), Y2.str()); printf ("actually calculating yields:"); if (setjmp (ovfl_buf)) { BadCond (Serious, "trap on underflow.\n"); } else { V9 = POW (HInvrse, Y2); printf (" %s .\n", V9.str()); if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) { BadCond (Serious, "this is not between 0 and underflow\n"); printf (" threshold = %s .\n", UfThold.str()); } else if (!(V9 > UfThold * (One + E9))) printf ("This computed value is O.K.\n"); else { BadCond (Defect, "this is not between 0 and underflow\n"); printf (" threshold = %s .\n", UfThold.str()); } } /*=============================================*/ Milestone = 160; /*=============================================*/ Pause (); printf ("Searching for Overflow threshold:\n"); printf ("This may generate an error.\n"); Y = -CInvrse; V9 = HInvrse * Y; if (setjmp (ovfl_buf)) { I = 0; V9 = Y; goto overflow; } do { V = Y; Y = V9; V9 = HInvrse * Y; } while (V9 < Y); I = 1; overflow: Z = V9; printf ("Can `Z = -Y' overflow?\n"); printf ("Trying it on Y = %s .\n", Y.str()); V9 = -Y; V0 = V9; if (V - Y == V + V0) printf ("Seems O.K.\n"); else { printf ("finds a "); BadCond (Flaw, "-(-Y) differs from Y.\n"); } if (Z != Y) { BadCond (Serious, ""); printf ("overflow past %s\n\tshrinks to %s .\n", Y.str(), Z.str()); } if (I) { Y = V * (HInvrse * U2 - HInvrse); Z = Y + ((One - HInvrse) * U2) * V; if (Z < V0) Y = Z; if (Y < V0) V = Y; if (V0 - V < V0) V = V0; } else { V = Y * (HInvrse * U2 - HInvrse); V = V + ((One - HInvrse) * U2) * Y; } printf ("Overflow threshold is V = %s .\n", V.str()); if (I) printf ("Overflow saturates at V0 = %s .\n", V0.str()); else printf ("There is no saturation value because " "the system traps on overflow.\n"); V9 = V * One; printf ("No Overflow should be signaled for V * 1 = %s\n", V9.str()); V9 = V / One; printf (" nor for V / 1 = %s.\n", V9.str()); printf ("Any overflow signal separating this * from the one\n"); printf ("above is a DEFECT.\n"); /*=============================================*/ Milestone = 170; /*=============================================*/ if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) { BadCond (Failure, "Comparisons involving "); printf ("+-%s, +-%s\nand +-%s are confused by Overflow.", V.str(), V0.str(), UfThold.str()); } /*=============================================*/ Milestone = 175; /*=============================================*/ printf ("\n"); for (Indx = 1; Indx <= 3; ++Indx) { switch (Indx) { case 1: Z = UfThold; break; case 2: Z = E0; break; case 3: Z = PseudoZero; break; } if (Z != Zero) { V9 = SQRT (Z); Y = V9 * V9; if (Y / (One - Radix * E9) < Z || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */ if (V9 > U1) BadCond (Serious, ""); else BadCond (Defect, ""); printf ("Comparison alleges that what prints as Z = %s\n", Z.str()); printf (" is too far from sqrt(Z) ^ 2 = %s .\n", Y.str()); } } } /*=============================================*/ Milestone = 180; /*=============================================*/ for (Indx = 1; Indx <= 2; ++Indx) { if (Indx == 1) Z = V; else Z = V0; V9 = SQRT (Z); X = (One - Radix * E9) * V9; V9 = V9 * X; if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) { Y = V9; if (X < W) BadCond (Serious, ""); else BadCond (Defect, ""); printf ("Comparison alleges that Z = %s\n", Z.str()); printf (" is too far from sqrt(Z) ^ 2 (%s) .\n", Y.str()); } } /*=============================================*/ Milestone = 190; /*=============================================*/ Pause (); X = UfThold * V; Y = Radix * Radix; if (X * Y < One || X > Y) { if (X * Y < U1 || X > Y / U1) BadCond (Defect, "Badly"); else BadCond (Flaw, ""); printf (" unbalanced range; UfThold * V = %s\n\t%s\n", X.str(), "is too far from 1.\n"); } /*=============================================*/ Milestone = 200; /*=============================================*/ for (Indx = 1; Indx <= 5; ++Indx) { X = F9; switch (Indx) { case 2: X = One + U2; break; case 3: X = V; break; case 4: X = UfThold; break; case 5: X = Radix; } Y = X; if (setjmp (ovfl_buf)) printf (" X / X traps when X = %s\n", X.str()); else { V9 = (Y / X - Half) - Half; if (V9 == Zero) continue; if (V9 == -U1 && Indx < 5) BadCond (Flaw, ""); else BadCond (Serious, ""); printf (" X / X differs from 1 when X = %s\n", X.str()); printf (" instead, X / X - 1/2 - 1/2 = %s .\n", V9.str()); } } /*=============================================*/ Milestone = 210; /*=============================================*/ MyZero = Zero; printf ("\n"); printf ("What message and/or values does Division by Zero produce?\n"); printf (" Trying to compute 1 / 0 produces ..."); if (!setjmp (ovfl_buf)) printf (" %s .\n", (One / MyZero).str()); printf ("\n Trying to compute 0 / 0 produces ..."); if (!setjmp (ovfl_buf)) printf (" %s .\n", (Zero / MyZero).str()); /*=============================================*/ Milestone = 220; /*=============================================*/ Pause (); printf ("\n"); { static const 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 ("Potentially 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 { if (StickyBit >= One && (Radix - Two) * (Radix - Nine - One) == Zero) { printf ("Rounding appears to conform to "); printf ("the proposed IEEE standard P"); if ((Radix == Two) && ((Precision - Four * Three * Two) * (Precision - TwentySeven - TwentySeven + One) == Zero)) 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"); } } printf ("END OF TEST.\n"); return 0; } template<typename FLOAT> FLOAT Paranoia<FLOAT>::Sign (FLOAT X) { return X >= FLOAT (long (0)) ? 1 : -1; } template<typename FLOAT> void Paranoia<FLOAT>::Pause () { if (do_pause) { fputs ("Press return...", stdout); fflush (stdout); getchar(); } printf ("\nDiagnosis resumes after milestone Number %d", Milestone); printf (" Page: %d\n\n", PageNo); ++Milestone; ++PageNo; } template<typename FLOAT> void Paranoia<FLOAT>::TstCond (int K, int Valid, const char *T) { if (!Valid) { BadCond (K, T); printf (".\n"); } } template<typename FLOAT> void Paranoia<FLOAT>::BadCond (int K, const char *T) { static const char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" }; ErrCnt[K] = ErrCnt[K] + 1; printf ("%s: %s", msg[K], T); } /* Random computes X = (Random1 + Random9)^5 Random1 = X - FLOOR(X) + 0.000005 * X; and returns the new value of Random1. */ template<typename FLOAT> FLOAT Paranoia<FLOAT>::Random () { FLOAT X, Y; X = Random1 + Random9; Y = X * X; Y = Y * Y; X = X * Y; Y = X - FLOOR (X); Random1 = Y + X * FLOAT ("0.000005"); return (Random1); } template<typename FLOAT> void Paranoia<FLOAT>::SqXMinX (int ErrKind) { FLOAT XA, XB; XB = X * BInvrse; XA = X - XB; SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp; if (SqEr != Zero) { if (SqEr < MinSqEr) MinSqEr = SqEr; if (SqEr > MaxSqEr) MaxSqEr = SqEr; J = J + 1; BadCond (ErrKind, "\n"); printf ("sqrt(%s) - %s = %s\n", (X * X).str(), X.str(), (OneUlp * SqEr).str()); printf ("\tinstead of correct value 0 .\n"); } } template<typename FLOAT> void Paranoia<FLOAT>::NewD () { X = Z1 * Q; X = FLOOR (Half - X / Radix) * Radix + X; Q = (Q - X * Z) / Radix + X * X * (D / Radix); Z = Z - Two * X * D; if (Z <= Zero) { Z = -Z; Z1 = -Z1; } D = Radix * D; } template<typename FLOAT> void Paranoia<FLOAT>::SR3750 () { if (!((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) { I = I + 1; X2 = SQRT (X * D); Y2 = (X2 - Z2) - (Y - Z2); X2 = X8 / (Y - Half); X2 = X2 - Half * X2 * X2; SqEr = (Y2 + Half) + (Half - X2); if (SqEr < MinSqEr) MinSqEr = SqEr; SqEr = Y2 - X2; if (SqEr > MaxSqEr) MaxSqEr = SqEr; } } template<typename FLOAT> void Paranoia<FLOAT>::IsYeqX () { if (Y != X) { if (N <= 0) { if (Z == Zero && Q <= Zero) printf ("WARNING: computing\n"); else BadCond (Defect, "computing\n"); printf ("\t(%s) ^ (%s)\n", Z.str(), Q.str()); printf ("\tyielded %s;\n", Y.str()); printf ("\twhich compared unequal to correct %s ;\n", X.str()); printf ("\t\tthey differ by %s .\n", (Y - X).str()); } N = N + 1; /* ... count discrepancies. */ } } template<typename FLOAT> void Paranoia<FLOAT>::PrintIfNPositive () { if (N > 0) printf ("Similar discrepancies have occurred %d times.\n", N); } template<typename FLOAT> void Paranoia<FLOAT>::TstPtUf () { N = 0; if (Z != Zero) { printf ("Since comparison denies Z = 0, evaluating "); printf ("(Z + Z) / Z should be safe.\n"); if (setjmp (ovfl_buf)) goto very_serious; Q9 = (Z + Z) / Z; printf ("What the machine gets for (Z + Z) / Z is %s .\n", Q9.str()); if (FABS (Q9 - Two) < Radix * U2) { printf ("This is O.K., provided Over/Underflow"); printf (" has NOT just been signaled.\n"); } else { if ((Q9 < One) || (Q9 > Two)) { very_serious: N = 1; ErrCnt[Serious] = ErrCnt[Serious] + 1; printf ("This is a VERY SERIOUS DEFECT!\n"); } else { N = 1; ErrCnt[Defect] = ErrCnt[Defect] + 1; printf ("This is a DEFECT!\n"); } } V9 = Z * One; Random1 = V9; V9 = One * Z; Random2 = V9; V9 = Z / One; if ((Z == Random1) && (Z == Random2) && (Z == V9)) { if (N > 0) Pause (); } else { N = 1; BadCond (Defect, "What prints as Z = "); printf ("%s\n\tcompares different from ", Z.str()); if (Z != Random1) printf ("Z * 1 = %s ", Random1.str()); if (!((Z == Random2) || (Random2 == Random1))) printf ("1 * Z == %s\n", Random2.str()); if (!(Z == V9)) printf ("Z / 1 = %s\n", V9.str()); if (Random2 != Random1) { ErrCnt[Defect] = ErrCnt[Defect] + 1; BadCond (Defect, "Multiplication does not commute!\n"); printf ("\tComparison alleges that 1 * Z = %s\n", Random2.str()); printf ("\tdiffers from Z * 1 = %s\n", Random1.str()); } Pause (); } } } template<typename FLOAT> void Paranoia<FLOAT>::notify (const char *s) { printf ("%s test appears to be inconsistent...\n", s); printf (" PLEASE NOTIFY KARPINKSI!\n"); } /* ====================================================================== */ int main(int ac, char **av) { setbuf(stdout, NULL); setbuf(stderr, NULL); while (1) switch (getopt (ac, av, "pvg:fdl")) { case -1: return 0; case 'p': do_pause = true; break; case 'v': verbose = true; break; case 'g': { static const struct { const char *name; const struct real_format *fmt; } fmts[] = { #define F(x) { #x, &x##_format } F(ieee_single), F(ieee_double), F(ieee_extended_motorola), F(ieee_extended_intel_96), F(ieee_extended_intel_128), F(ibm_extended), F(ieee_quad), F(vax_f), F(vax_d), F(vax_g), F(i370_single), F(i370_double), F(c4x_single), F(c4x_extended), F(real_internal), #undef F }; int i, n = sizeof (fmts)/sizeof(*fmts); for (i = 0; i < n; ++i) if (strcmp (fmts[i].name, optarg) == 0) break; if (i == n) { printf ("Unknown implementation \"%s\"; " "available implementations:\n", optarg); for (i = 0; i < n; ++i) printf ("\t%s\n", fmts[i].name); return 1; } // We cheat and use the same mode all the time, but vary // the format used for that mode. real_format_for_mode[int(real_c_float::MODE) - int(QFmode)] = fmts[i].fmt; Paranoia<real_c_float>().main(); break; } case 'f': Paranoia < native_float<float> >().main(); break; case 'd': Paranoia < native_float<double> >().main(); break; case 'l': #ifndef NO_LONG_DOUBLE Paranoia < native_float<long double> >().main(); #endif break; case '?': puts ("-p\tpause between pages"); puts ("-g<FMT>\treal.c implementation FMT"); puts ("-f\tnative float"); puts ("-d\tnative double"); puts ("-l\tnative long double"); return 0; } } /* GCC stuff referenced by real.o. */ extern "C" void fancy_abort () { abort (); } int target_flags = 0; extern "C" int floor_log2_wide (unsigned HOST_WIDE_INT x) { int log = -1; while (x != 0) log++, x >>= 1; return log; }
Go to most recent revision | Compare with Previous | Blame | View Log