OpenE2K
/
gcc
2
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
gcc/contrib/paranoia.cc

2715 lines
65 KiB
C++
Raw Blame History

/* 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(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;
}
Reference in New Issue View Git Blame Copy Permalink