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