mirror of
git://gcc.gnu.org/git/gcc.git
synced 2024-12-21 15:19:46 +08:00
ad03e4c4c7
* paranoia.cc (real_c_float::image): Accomodate size of
real_internal format.
(main): Unbuffer stdio. Add real_internal format.
(floor_log2_wide): New.
From-SVN: r58395
2717 lines
65 KiB
C++
2717 lines
65 KiB
C++
/* 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 excersize. Templated so that we can
|
|
all fp sizes. */
|
|
|
|
class real_c_float
|
|
{
|
|
public:
|
|
static const enum machine_mode MODE = SFmode;
|
|
|
|
private:
|
|
static const int external_max = 128 / 32;
|
|
static const int internal_max
|
|
= (sizeof (REAL_VALUE_TYPE) + sizeof (long) + 1) / sizeof (long);
|
|
long image[external_max < internal_max ? internal_max : external_max];
|
|
|
|
void from_long(long);
|
|
void from_str(const char *);
|
|
void binop(int code, const real_c_float&);
|
|
void unop(int code);
|
|
bool cmp(int code, const real_c_float&) const;
|
|
|
|
public:
|
|
real_c_float()
|
|
{ }
|
|
real_c_float(long l)
|
|
{ from_long(l); }
|
|
real_c_float(const char *s)
|
|
{ from_str(s); }
|
|
real_c_float(const real_c_float &b)
|
|
{ memcpy(image, b.image, sizeof(image)); }
|
|
|
|
const real_c_float& operator= (long l)
|
|
{ from_long(l); return *this; }
|
|
const real_c_float& operator= (const char *s)
|
|
{ from_str(s); return *this; }
|
|
const real_c_float& operator= (const real_c_float &b)
|
|
{ memcpy(image, b.image, sizeof(image)); return *this; }
|
|
|
|
const real_c_float& operator+= (const real_c_float &b)
|
|
{ binop(PLUS_EXPR, b); return *this; }
|
|
const real_c_float& operator-= (const real_c_float &b)
|
|
{ binop(MINUS_EXPR, b); return *this; }
|
|
const real_c_float& operator*= (const real_c_float &b)
|
|
{ binop(MULT_EXPR, b); return *this; }
|
|
const real_c_float& operator/= (const real_c_float &b)
|
|
{ binop(RDIV_EXPR, b); return *this; }
|
|
|
|
real_c_float operator- () const
|
|
{ real_c_float r(*this); r.unop(NEGATE_EXPR); return r; }
|
|
real_c_float abs () const
|
|
{ real_c_float r(*this); r.unop(ABS_EXPR); return r; }
|
|
|
|
bool operator < (const real_c_float &b) const { return cmp(LT_EXPR, b); }
|
|
bool operator <= (const real_c_float &b) const { return cmp(LE_EXPR, b); }
|
|
bool operator == (const real_c_float &b) const { return cmp(EQ_EXPR, b); }
|
|
bool operator != (const real_c_float &b) const { return cmp(NE_EXPR, b); }
|
|
bool operator >= (const real_c_float &b) const { return cmp(GE_EXPR, b); }
|
|
bool operator > (const real_c_float &b) const { return cmp(GT_EXPR, b); }
|
|
|
|
const char * str () const;
|
|
const char * hex () const;
|
|
long integer () const;
|
|
int exp () const;
|
|
void ldexp (int);
|
|
};
|
|
|
|
void
|
|
real_c_float::from_long (long l)
|
|
{
|
|
REAL_VALUE_TYPE f;
|
|
|
|
real_from_integer (&f, MODE, l, l < 0 ? -1 : 0, 0);
|
|
real_to_target (image, &f, MODE);
|
|
}
|
|
|
|
void
|
|
real_c_float::from_str (const char *s)
|
|
{
|
|
REAL_VALUE_TYPE f;
|
|
const char *p = s;
|
|
|
|
if (*p == '-' || *p == '+')
|
|
p++;
|
|
if (strcasecmp(p, "inf") == 0)
|
|
{
|
|
real_inf (&f);
|
|
if (*s == '-')
|
|
real_arithmetic (&f, NEGATE_EXPR, &f, NULL);
|
|
}
|
|
else if (strcasecmp(p, "nan") == 0)
|
|
real_nan (&f, "", 1, MODE);
|
|
else
|
|
real_from_string (&f, s);
|
|
|
|
real_to_target (image, &f, MODE);
|
|
}
|
|
|
|
void
|
|
real_c_float::binop (int code, const real_c_float &b)
|
|
{
|
|
REAL_VALUE_TYPE ai, bi, ri;
|
|
|
|
real_from_target (&ai, image, MODE);
|
|
real_from_target (&bi, b.image, MODE);
|
|
real_arithmetic (&ri, code, &ai, &bi);
|
|
real_to_target (image, &ri, MODE);
|
|
|
|
if (verbose)
|
|
{
|
|
char ab[64], bb[64], rb[64];
|
|
const real_format *fmt = real_format_for_mode[MODE - QFmode];
|
|
const int digits = (fmt->p * fmt->log2_b + 3) / 4;
|
|
char symbol_for_code;
|
|
|
|
real_from_target (&ri, image, MODE);
|
|
real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
|
|
real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0);
|
|
real_to_hexadecimal (rb, &ri, sizeof(rb), digits, 0);
|
|
|
|
switch (code)
|
|
{
|
|
case PLUS_EXPR:
|
|
symbol_for_code = '+';
|
|
break;
|
|
case MINUS_EXPR:
|
|
symbol_for_code = '-';
|
|
break;
|
|
case MULT_EXPR:
|
|
symbol_for_code = '*';
|
|
break;
|
|
case RDIV_EXPR:
|
|
symbol_for_code = '/';
|
|
break;
|
|
default:
|
|
abort ();
|
|
}
|
|
|
|
fprintf (stderr, "%6d: %s %c %s = %s\n", verbose_index++,
|
|
ab, symbol_for_code, bb, rb);
|
|
}
|
|
}
|
|
|
|
void
|
|
real_c_float::unop (int code)
|
|
{
|
|
REAL_VALUE_TYPE ai, ri;
|
|
|
|
real_from_target (&ai, image, MODE);
|
|
real_arithmetic (&ri, code, &ai, NULL);
|
|
real_to_target (image, &ri, MODE);
|
|
|
|
if (verbose)
|
|
{
|
|
char ab[64], rb[64];
|
|
const real_format *fmt = real_format_for_mode[MODE - QFmode];
|
|
const int digits = (fmt->p * fmt->log2_b + 3) / 4;
|
|
const char *symbol_for_code;
|
|
|
|
real_from_target (&ri, image, MODE);
|
|
real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
|
|
real_to_hexadecimal (rb, &ri, sizeof(rb), digits, 0);
|
|
|
|
switch (code)
|
|
{
|
|
case NEGATE_EXPR:
|
|
symbol_for_code = "-";
|
|
break;
|
|
case ABS_EXPR:
|
|
symbol_for_code = "abs ";
|
|
break;
|
|
default:
|
|
abort ();
|
|
}
|
|
|
|
fprintf (stderr, "%6d: %s%s = %s\n", verbose_index++,
|
|
symbol_for_code, ab, rb);
|
|
}
|
|
}
|
|
|
|
bool
|
|
real_c_float::cmp (int code, const real_c_float &b) const
|
|
{
|
|
REAL_VALUE_TYPE ai, bi;
|
|
bool ret;
|
|
|
|
real_from_target (&ai, image, MODE);
|
|
real_from_target (&bi, b.image, MODE);
|
|
ret = real_compare (code, &ai, &bi);
|
|
|
|
if (verbose)
|
|
{
|
|
char ab[64], bb[64];
|
|
const real_format *fmt = real_format_for_mode[MODE - QFmode];
|
|
const int digits = (fmt->p * fmt->log2_b + 3) / 4;
|
|
const char *symbol_for_code;
|
|
|
|
real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
|
|
real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0);
|
|
|
|
switch (code)
|
|
{
|
|
case LT_EXPR:
|
|
symbol_for_code = "<";
|
|
break;
|
|
case LE_EXPR:
|
|
symbol_for_code = "<=";
|
|
break;
|
|
case EQ_EXPR:
|
|
symbol_for_code = "==";
|
|
break;
|
|
case NE_EXPR:
|
|
symbol_for_code = "!=";
|
|
break;
|
|
case GE_EXPR:
|
|
symbol_for_code = ">=";
|
|
break;
|
|
case GT_EXPR:
|
|
symbol_for_code = ">";
|
|
break;
|
|
default:
|
|
abort ();
|
|
}
|
|
|
|
fprintf (stderr, "%6d: %s %s %s = %s\n", verbose_index++,
|
|
ab, symbol_for_code, bb, (ret ? "true" : "false"));
|
|
}
|
|
|
|
return ret;
|
|
}
|
|
|
|
const char *
|
|
real_c_float::str() const
|
|
{
|
|
REAL_VALUE_TYPE f;
|
|
const real_format *fmt = real_format_for_mode[MODE - QFmode];
|
|
const int digits = int(fmt->p * fmt->log2_b * .30102999566398119521 + 1);
|
|
|
|
real_from_target (&f, image, MODE);
|
|
char *buf = new char[digits + 10];
|
|
real_to_decimal (buf, &f, digits+10, digits, 0);
|
|
|
|
return buf;
|
|
}
|
|
|
|
const char *
|
|
real_c_float::hex() const
|
|
{
|
|
REAL_VALUE_TYPE f;
|
|
const real_format *fmt = real_format_for_mode[MODE - QFmode];
|
|
const int digits = (fmt->p * fmt->log2_b + 3) / 4;
|
|
|
|
real_from_target (&f, image, MODE);
|
|
char *buf = new char[digits + 10];
|
|
real_to_hexadecimal (buf, &f, digits+10, digits, 0);
|
|
|
|
return buf;
|
|
}
|
|
|
|
long
|
|
real_c_float::integer() const
|
|
{
|
|
REAL_VALUE_TYPE f;
|
|
real_from_target (&f, image, MODE);
|
|
return real_to_integer (&f);
|
|
}
|
|
|
|
int
|
|
real_c_float::exp() const
|
|
{
|
|
REAL_VALUE_TYPE f;
|
|
real_from_target (&f, image, MODE);
|
|
return real_exponent (&f);
|
|
}
|
|
|
|
void
|
|
real_c_float::ldexp (int exp)
|
|
{
|
|
REAL_VALUE_TYPE ai;
|
|
|
|
real_from_target (&ai, image, MODE);
|
|
real_ldexp (&ai, &ai, exp);
|
|
real_to_target (image, &ai, MODE);
|
|
}
|
|
|
|
/* ====================================================================== */
|
|
/* An implementation of the abstract floating point class that uses native
|
|
arithmetic. Exists for reference and debugging. */
|
|
|
|
template<typename T>
|
|
class native_float
|
|
{
|
|
private:
|
|
// Force intermediate results back to memory.
|
|
volatile T image;
|
|
|
|
static T from_str (const char *);
|
|
static T do_abs (T);
|
|
static T verbose_binop (T, char, T, T);
|
|
static T verbose_unop (const char *, T, T);
|
|
static bool verbose_cmp (T, const char *, T, bool);
|
|
|
|
public:
|
|
native_float()
|
|
{ }
|
|
native_float(long l)
|
|
{ image = l; }
|
|
native_float(const char *s)
|
|
{ image = from_str(s); }
|
|
native_float(const native_float &b)
|
|
{ image = b.image; }
|
|
|
|
const native_float& operator= (long l)
|
|
{ image = l; return *this; }
|
|
const native_float& operator= (const char *s)
|
|
{ image = from_str(s); return *this; }
|
|
const native_float& operator= (const native_float &b)
|
|
{ image = b.image; return *this; }
|
|
|
|
const native_float& operator+= (const native_float &b)
|
|
{
|
|
image = verbose_binop(image, '+', b.image, image + b.image);
|
|
return *this;
|
|
}
|
|
const native_float& operator-= (const native_float &b)
|
|
{
|
|
image = verbose_binop(image, '-', b.image, image - b.image);
|
|
return *this;
|
|
}
|
|
const native_float& operator*= (const native_float &b)
|
|
{
|
|
image = verbose_binop(image, '*', b.image, image * b.image);
|
|
return *this;
|
|
}
|
|
const native_float& operator/= (const native_float &b)
|
|
{
|
|
image = verbose_binop(image, '/', b.image, image / b.image);
|
|
return *this;
|
|
}
|
|
|
|
native_float operator- () const
|
|
{
|
|
native_float r;
|
|
r.image = verbose_unop("-", image, -image);
|
|
return r;
|
|
}
|
|
native_float abs () const
|
|
{
|
|
native_float r;
|
|
r.image = verbose_unop("abs ", image, do_abs(image));
|
|
return r;
|
|
}
|
|
|
|
bool operator < (const native_float &b) const
|
|
{ return verbose_cmp(image, "<", b.image, image < b.image); }
|
|
bool operator <= (const native_float &b) const
|
|
{ return verbose_cmp(image, "<=", b.image, image <= b.image); }
|
|
bool operator == (const native_float &b) const
|
|
{ return verbose_cmp(image, "==", b.image, image == b.image); }
|
|
bool operator != (const native_float &b) const
|
|
{ return verbose_cmp(image, "!=", b.image, image != b.image); }
|
|
bool operator >= (const native_float &b) const
|
|
{ return verbose_cmp(image, ">=", b.image, image >= b.image); }
|
|
bool operator > (const native_float &b) const
|
|
{ return verbose_cmp(image, ">", b.image, image > b.image); }
|
|
|
|
const char * str () const;
|
|
const char * hex () const;
|
|
long integer () const
|
|
{ return long(image); }
|
|
int exp () const;
|
|
void ldexp (int);
|
|
};
|
|
|
|
template<typename T>
|
|
inline T
|
|
native_float<T>::from_str (const char *s)
|
|
{
|
|
return strtold (s, NULL);
|
|
}
|
|
|
|
template<>
|
|
inline float
|
|
native_float<float>::from_str (const char *s)
|
|
{
|
|
return strtof (s, NULL);
|
|
}
|
|
|
|
template<>
|
|
inline double
|
|
native_float<double>::from_str (const char *s)
|
|
{
|
|
return strtod (s, NULL);
|
|
}
|
|
|
|
template<typename T>
|
|
inline T
|
|
native_float<T>::do_abs (T image)
|
|
{
|
|
return fabsl (image);
|
|
}
|
|
|
|
template<>
|
|
inline float
|
|
native_float<float>::do_abs (float image)
|
|
{
|
|
return fabsf (image);
|
|
}
|
|
|
|
template<>
|
|
inline double
|
|
native_float<double>::do_abs (double image)
|
|
{
|
|
return fabs (image);
|
|
}
|
|
|
|
template<typename T>
|
|
T
|
|
native_float<T>::verbose_binop (T a, char symbol, T b, T r)
|
|
{
|
|
if (verbose)
|
|
{
|
|
const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1;
|
|
#ifdef NO_LONG_DOUBLE
|
|
fprintf (stderr, "%6d: %.*a %c %.*a = %.*a\n", verbose_index++,
|
|
digits, (double)a, symbol,
|
|
digits, (double)b, digits, (double)r);
|
|
#else
|
|
fprintf (stderr, "%6d: %.*La %c %.*La = %.*La\n", verbose_index++,
|
|
digits, (long double)a, symbol,
|
|
digits, (long double)b, digits, (long double)r);
|
|
#endif
|
|
}
|
|
return r;
|
|
}
|
|
|
|
template<typename T>
|
|
T
|
|
native_float<T>::verbose_unop (const char *symbol, T a, T r)
|
|
{
|
|
if (verbose)
|
|
{
|
|
const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1;
|
|
#ifdef NO_LONG_DOUBLE
|
|
fprintf (stderr, "%6d: %s%.*a = %.*a\n", verbose_index++,
|
|
symbol, digits, (double)a, digits, (double)r);
|
|
#else
|
|
fprintf (stderr, "%6d: %s%.*La = %.*La\n", verbose_index++,
|
|
symbol, digits, (long double)a, digits, (long double)r);
|
|
#endif
|
|
}
|
|
return r;
|
|
}
|
|
|
|
template<typename T>
|
|
bool
|
|
native_float<T>::verbose_cmp (T a, const char *symbol, T b, bool r)
|
|
{
|
|
if (verbose)
|
|
{
|
|
const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1;
|
|
#ifndef NO_LONG_DOUBLE
|
|
fprintf (stderr, "%6d: %.*a %s %.*a = %s\n", verbose_index++,
|
|
digits, (double)a, symbol,
|
|
digits, (double)b, (r ? "true" : "false"));
|
|
#else
|
|
fprintf (stderr, "%6d: %.*La %s %.*La = %s\n", verbose_index++,
|
|
digits, (long double)a, symbol,
|
|
digits, (long double)b, (r ? "true" : "false"));
|
|
#endif
|
|
}
|
|
return r;
|
|
}
|
|
|
|
template<typename T>
|
|
const char *
|
|
native_float<T>::str() const
|
|
{
|
|
char *buf = new char[50];
|
|
const int digits = int(sizeof(T) * CHAR_BIT * .30102999566398119521 + 1);
|
|
#ifndef NO_LONG_DOUBLE
|
|
sprintf (buf, "%.*e", digits - 1, (double) image);
|
|
#else
|
|
sprintf (buf, "%.*Le", digits - 1, (long double) image);
|
|
#endif
|
|
return buf;
|
|
}
|
|
|
|
template<typename T>
|
|
const char *
|
|
native_float<T>::hex() const
|
|
{
|
|
char *buf = new char[50];
|
|
const int digits = int(sizeof(T) * CHAR_BIT / 4);
|
|
#ifndef NO_LONG_DOUBLE
|
|
sprintf (buf, "%.*a", digits - 1, (double) image);
|
|
#else
|
|
sprintf (buf, "%.*La", digits - 1, (long double) image);
|
|
#endif
|
|
return buf;
|
|
}
|
|
|
|
template<typename T>
|
|
int
|
|
native_float<T>::exp() const
|
|
{
|
|
int e;
|
|
frexp (image, &e);
|
|
return e;
|
|
}
|
|
|
|
template<typename T>
|
|
void
|
|
native_float<T>::ldexp (int exp)
|
|
{
|
|
image = ldexpl (image, exp);
|
|
}
|
|
|
|
template<>
|
|
void
|
|
native_float<float>::ldexp (int exp)
|
|
{
|
|
image = ldexpf (image, exp);
|
|
}
|
|
|
|
template<>
|
|
void
|
|
native_float<double>::ldexp (int exp)
|
|
{
|
|
image = ::ldexp (image, exp);
|
|
}
|
|
|
|
/* ====================================================================== */
|
|
/* Some libm routines that Paranoia expects to be available. */
|
|
|
|
template<typename FLOAT>
|
|
inline FLOAT
|
|
FABS (const FLOAT &f)
|
|
{
|
|
return f.abs();
|
|
}
|
|
|
|
template<typename FLOAT, typename RHS>
|
|
inline FLOAT
|
|
operator+ (const FLOAT &a, const RHS &b)
|
|
{
|
|
return FLOAT(a) += FLOAT(b);
|
|
}
|
|
|
|
template<typename FLOAT, typename RHS>
|
|
inline FLOAT
|
|
operator- (const FLOAT &a, const RHS &b)
|
|
{
|
|
return FLOAT(a) -= FLOAT(b);
|
|
}
|
|
|
|
template<typename FLOAT, typename RHS>
|
|
inline FLOAT
|
|
operator* (const FLOAT &a, const RHS &b)
|
|
{
|
|
return FLOAT(a) *= FLOAT(b);
|
|
}
|
|
|
|
template<typename FLOAT, typename RHS>
|
|
inline FLOAT
|
|
operator/ (const FLOAT &a, const RHS &b)
|
|
{
|
|
return FLOAT(a) /= FLOAT(b);
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
FLOOR (const FLOAT &f)
|
|
{
|
|
/* ??? This is only correct when F is representable as an integer. */
|
|
long i = f.integer();
|
|
FLOAT r;
|
|
|
|
r = i;
|
|
if (i < 0 && f != r)
|
|
r = i - 1;
|
|
|
|
return r;
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
SQRT (const FLOAT &f)
|
|
{
|
|
#if 0
|
|
FLOAT zero = long(0);
|
|
FLOAT two = 2;
|
|
FLOAT one = 1;
|
|
FLOAT diff, diff2;
|
|
FLOAT z, t;
|
|
|
|
if (f == zero)
|
|
return zero;
|
|
if (f < zero)
|
|
return zero / zero;
|
|
if (f == one)
|
|
return f;
|
|
|
|
z = f;
|
|
z.ldexp (-f.exp() / 2);
|
|
|
|
diff2 = FABS (z * z - f);
|
|
if (diff2 > zero)
|
|
while (1)
|
|
{
|
|
t = (f / (two * z)) + (z / two);
|
|
diff = FABS (t * t - f);
|
|
if (diff >= diff2)
|
|
break;
|
|
z = t;
|
|
diff2 = diff;
|
|
}
|
|
|
|
return z;
|
|
#elif defined(NO_LONG_DOUBLE)
|
|
double d;
|
|
char buf[64];
|
|
|
|
d = strtod (f.hex(), NULL);
|
|
d = sqrt (d);
|
|
sprintf(buf, "%.35a", d);
|
|
|
|
return FLOAT(buf);
|
|
#else
|
|
long double ld;
|
|
char buf[64];
|
|
|
|
ld = strtold (f.hex(), NULL);
|
|
ld = sqrtl (ld);
|
|
sprintf(buf, "%.35La", ld);
|
|
|
|
return FLOAT(buf);
|
|
#endif
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
LOG (FLOAT x)
|
|
{
|
|
#if 0
|
|
FLOAT zero = long(0);
|
|
FLOAT one = 1;
|
|
|
|
if (x <= zero)
|
|
return zero / zero;
|
|
if (x == one)
|
|
return zero;
|
|
|
|
int exp = x.exp() - 1;
|
|
x.ldexp(-exp);
|
|
|
|
FLOAT xm1 = x - one;
|
|
FLOAT y = xm1;
|
|
long n = 2;
|
|
|
|
FLOAT sum = xm1;
|
|
while (1)
|
|
{
|
|
y *= xm1;
|
|
FLOAT term = y / FLOAT (n);
|
|
FLOAT next = sum + term;
|
|
if (next == sum)
|
|
break;
|
|
sum = next;
|
|
if (++n == 1000)
|
|
break;
|
|
}
|
|
|
|
if (exp)
|
|
sum += FLOAT (exp) * FLOAT(".69314718055994530941");
|
|
|
|
return sum;
|
|
#elif defined (NO_LONG_DOUBLE)
|
|
double d;
|
|
char buf[64];
|
|
|
|
d = strtod (x.hex(), NULL);
|
|
d = log (d);
|
|
sprintf(buf, "%.35a", d);
|
|
|
|
return FLOAT(buf);
|
|
#else
|
|
long double ld;
|
|
char buf[64];
|
|
|
|
ld = strtold (x.hex(), NULL);
|
|
ld = logl (ld);
|
|
sprintf(buf, "%.35La", ld);
|
|
|
|
return FLOAT(buf);
|
|
#endif
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
EXP (const FLOAT &x)
|
|
{
|
|
/* Cheat. */
|
|
#ifdef NO_LONG_DOUBLE
|
|
double d;
|
|
char buf[64];
|
|
|
|
d = strtod (x.hex(), NULL);
|
|
d = exp (d);
|
|
sprintf(buf, "%.35a", d);
|
|
|
|
return FLOAT(buf);
|
|
#else
|
|
long double ld;
|
|
char buf[64];
|
|
|
|
ld = strtold (x.hex(), NULL);
|
|
ld = expl (ld);
|
|
sprintf(buf, "%.35La", ld);
|
|
|
|
return FLOAT(buf);
|
|
#endif
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
POW (const FLOAT &base, const FLOAT &exp)
|
|
{
|
|
/* Cheat. */
|
|
#ifdef NO_LONG_DOUBLE
|
|
double d1, d2;
|
|
char buf[64];
|
|
|
|
d1 = strtod (base.hex(), NULL);
|
|
d2 = strtod (exp.hex(), NULL);
|
|
d1 = pow (d1, d2);
|
|
sprintf(buf, "%.35a", d1);
|
|
|
|
return FLOAT(buf);
|
|
#else
|
|
long double ld1, ld2;
|
|
char buf[64];
|
|
|
|
ld1 = strtold (base.hex(), NULL);
|
|
ld2 = strtold (exp.hex(), NULL);
|
|
ld1 = powl (ld1, ld2);
|
|
sprintf(buf, "%.35La", ld1);
|
|
|
|
return FLOAT(buf);
|
|
#endif
|
|
}
|
|
|
|
/* ====================================================================== */
|
|
/* Real Paranoia begins again here. We wrap the thing in a template so
|
|
that we can instantiate it for each floating point type we care for. */
|
|
|
|
int NoTrials = 20; /*Number of tests for commutativity. */
|
|
bool do_pause = false;
|
|
|
|
enum Guard { No, Yes };
|
|
enum Rounding { Other, Rounded, Chopped };
|
|
enum Class { Failure, Serious, Defect, Flaw };
|
|
|
|
template<typename FLOAT>
|
|
struct Paranoia
|
|
{
|
|
FLOAT Radix, BInvrse, RadixD2, BMinusU2;
|
|
|
|
/* Small floating point constants. */
|
|
FLOAT Zero;
|
|
FLOAT Half;
|
|
FLOAT One;
|
|
FLOAT Two;
|
|
FLOAT Three;
|
|
FLOAT Four;
|
|
FLOAT Five;
|
|
FLOAT Eight;
|
|
FLOAT Nine;
|
|
FLOAT TwentySeven;
|
|
FLOAT ThirtyTwo;
|
|
FLOAT TwoForty;
|
|
FLOAT MinusOne;
|
|
FLOAT OneAndHalf;
|
|
|
|
/* Declarations of Variables. */
|
|
int Indx;
|
|
char ch[8];
|
|
FLOAT AInvrse, A1;
|
|
FLOAT C, CInvrse;
|
|
FLOAT D, FourD;
|
|
FLOAT E0, E1, Exp2, E3, MinSqEr;
|
|
FLOAT SqEr, MaxSqEr, E9;
|
|
FLOAT Third;
|
|
FLOAT F6, F9;
|
|
FLOAT H, HInvrse;
|
|
int I;
|
|
FLOAT StickyBit, J;
|
|
FLOAT MyZero;
|
|
FLOAT Precision;
|
|
FLOAT Q, Q9;
|
|
FLOAT R, Random9;
|
|
FLOAT T, Underflow, S;
|
|
FLOAT OneUlp, UfThold, U1, U2;
|
|
FLOAT V, V0, V9;
|
|
FLOAT W;
|
|
FLOAT X, X1, X2, X8, Random1;
|
|
FLOAT Y, Y1, Y2, Random2;
|
|
FLOAT Z, PseudoZero, Z1, Z2, Z9;
|
|
int ErrCnt[4];
|
|
int Milestone;
|
|
int PageNo;
|
|
int M, N, N1;
|
|
Guard GMult, GDiv, GAddSub;
|
|
Rounding RMult, RDiv, RAddSub, RSqrt;
|
|
int Break, Done, NotMonot, Monot, Anomaly, IEEE, SqRWrng, UfNGrad;
|
|
|
|
/* Computed constants. */
|
|
/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
|
|
/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
|
|
|
|
int main ();
|
|
|
|
FLOAT Sign (FLOAT);
|
|
FLOAT Random ();
|
|
void Pause ();
|
|
void BadCond (int, const char *);
|
|
void SqXMinX (int);
|
|
void TstCond (int, int, const char *);
|
|
void notify (const char *);
|
|
void IsYeqX ();
|
|
void NewD ();
|
|
void PrintIfNPositive ();
|
|
void SR3750 ();
|
|
void TstPtUf ();
|
|
|
|
// Pretend we're bss.
|
|
Paranoia() { memset(this, 0, sizeof (*this)); }
|
|
};
|
|
|
|
template<typename FLOAT>
|
|
int
|
|
Paranoia<FLOAT>::main()
|
|
{
|
|
/* First two assignments use integer right-hand sides. */
|
|
Zero = long(0);
|
|
One = long(1);
|
|
Two = long(2);
|
|
Three = long(3);
|
|
Four = long(4);
|
|
Five = long(5);
|
|
Eight = long(8);
|
|
Nine = long(9);
|
|
TwentySeven = long(27);
|
|
ThirtyTwo = long(32);
|
|
TwoForty = long(240);
|
|
MinusOne = long(-1);
|
|
Half = "0x1p-1";
|
|
OneAndHalf = "0x3p-1";
|
|
ErrCnt[Failure] = 0;
|
|
ErrCnt[Serious] = 0;
|
|
ErrCnt[Defect] = 0;
|
|
ErrCnt[Flaw] = 0;
|
|
PageNo = 1;
|
|
/*=============================================*/
|
|
Milestone = 7;
|
|
/*=============================================*/
|
|
printf ("Program is now RUNNING tests on small integers:\n");
|
|
|
|
TstCond (Failure, (Zero + Zero == Zero), "0+0 != 0");
|
|
TstCond (Failure, (One - One == Zero), "1-1 != 0");
|
|
TstCond (Failure, (One > Zero), "1 <= 0");
|
|
TstCond (Failure, (One + One == Two), "1+1 != 2");
|
|
|
|
Z = -Zero;
|
|
if (Z != Zero)
|
|
{
|
|
ErrCnt[Failure] = ErrCnt[Failure] + 1;
|
|
printf ("Comparison alleges that -0.0 is Non-zero!\n");
|
|
U2 = "0.001";
|
|
Radix = 1;
|
|
TstPtUf ();
|
|
}
|
|
|
|
TstCond (Failure, (Three == Two + One), "3 != 2+1");
|
|
TstCond (Failure, (Four == Three + One), "4 != 3+1");
|
|
TstCond (Failure, (Four + Two * (-Two) == Zero), "4 + 2*(-2) != 0");
|
|
TstCond (Failure, (Four - Three - One == Zero), "4-3-1 != 0");
|
|
|
|
TstCond (Failure, (MinusOne == (Zero - One)), "-1 != 0-1");
|
|
TstCond (Failure, (MinusOne + One == Zero), "-1+1 != 0");
|
|
TstCond (Failure, (One + MinusOne == Zero), "1+(-1) != 0");
|
|
TstCond (Failure, (MinusOne + FABS (One) == Zero), "-1+abs(1) != 0");
|
|
TstCond (Failure, (MinusOne + MinusOne * MinusOne == Zero),
|
|
"-1+(-1)*(-1) != 0");
|
|
|
|
TstCond (Failure, Half + MinusOne + Half == Zero, "1/2 + (-1) + 1/2 != 0");
|
|
|
|
/*=============================================*/
|
|
Milestone = 10;
|
|
/*=============================================*/
|
|
|
|
TstCond (Failure, (Nine == Three * Three), "9 != 3*3");
|
|
TstCond (Failure, (TwentySeven == Nine * Three), "27 != 9*3");
|
|
TstCond (Failure, (Eight == Four + Four), "8 != 4+4");
|
|
TstCond (Failure, (ThirtyTwo == Eight * Four), "32 != 8*4");
|
|
TstCond (Failure, (ThirtyTwo - TwentySeven - Four - One == Zero),
|
|
"32-27-4-1 != 0");
|
|
|
|
TstCond (Failure, Five == Four + One, "5 != 4+1");
|
|
TstCond (Failure, TwoForty == Four * Five * Three * Four, "240 != 4*5*3*4");
|
|
TstCond (Failure, TwoForty / Three - Four * Four * Five == Zero,
|
|
"240/3 - 4*4*5 != 0");
|
|
TstCond (Failure, TwoForty / Four - Five * Three * Four == Zero,
|
|
"240/4 - 5*3*4 != 0");
|
|
TstCond (Failure, TwoForty / Five - Four * Three * Four == Zero,
|
|
"240/5 - 4*3*4 != 0");
|
|
|
|
if (ErrCnt[Failure] == 0)
|
|
{
|
|
printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
|
|
printf ("\n");
|
|
}
|
|
printf ("Searching for Radix and Precision.\n");
|
|
W = One;
|
|
do
|
|
{
|
|
W = W + W;
|
|
Y = W + One;
|
|
Z = Y - W;
|
|
Y = Z - One;
|
|
}
|
|
while (MinusOne + FABS (Y) < Zero);
|
|
/*.. now W is just big enough that |((W+1)-W)-1| >= 1 ... */
|
|
Precision = Zero;
|
|
Y = One;
|
|
do
|
|
{
|
|
Radix = W + Y;
|
|
Y = Y + Y;
|
|
Radix = Radix - W;
|
|
}
|
|
while (Radix == Zero);
|
|
if (Radix < Two)
|
|
Radix = One;
|
|
printf ("Radix = %s .\n", Radix.str());
|
|
if (Radix != One)
|
|
{
|
|
W = One;
|
|
do
|
|
{
|
|
Precision = Precision + One;
|
|
W = W * Radix;
|
|
Y = W + One;
|
|
}
|
|
while ((Y - W) == One);
|
|
}
|
|
/*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
|
|
... */
|
|
U1 = One / W;
|
|
U2 = Radix * U1;
|
|
printf ("Closest relative separation found is U1 = %s .\n\n", U1.str());
|
|
printf ("Recalculating radix and precision\n ");
|
|
|
|
/*save old values */
|
|
E0 = Radix;
|
|
E1 = U1;
|
|
E9 = U2;
|
|
E3 = Precision;
|
|
|
|
X = Four / Three;
|
|
Third = X - One;
|
|
F6 = Half - Third;
|
|
X = F6 + F6;
|
|
X = FABS (X - Third);
|
|
if (X < U2)
|
|
X = U2;
|
|
|
|
/*... now X = (unknown no.) ulps of 1+... */
|
|
do
|
|
{
|
|
U2 = X;
|
|
Y = Half * U2 + ThirtyTwo * U2 * U2;
|
|
Y = One + Y;
|
|
X = Y - One;
|
|
}
|
|
while (!((U2 <= X) || (X <= Zero)));
|
|
|
|
/*... now U2 == 1 ulp of 1 + ... */
|
|
X = Two / Three;
|
|
F6 = X - Half;
|
|
Third = F6 + F6;
|
|
X = Third - Half;
|
|
X = FABS (X + F6);
|
|
if (X < U1)
|
|
X = U1;
|
|
|
|
/*... now X == (unknown no.) ulps of 1 -... */
|
|
do
|
|
{
|
|
U1 = X;
|
|
Y = Half * U1 + ThirtyTwo * U1 * U1;
|
|
Y = Half - Y;
|
|
X = Half + Y;
|
|
Y = Half - X;
|
|
X = Half + Y;
|
|
}
|
|
while (!((U1 <= X) || (X <= Zero)));
|
|
/*... now U1 == 1 ulp of 1 - ... */
|
|
if (U1 == E1)
|
|
printf ("confirms closest relative separation U1 .\n");
|
|
else
|
|
printf ("gets better closest relative separation U1 = %s .\n", U1.str());
|
|
W = One / U1;
|
|
F9 = (Half - U1) + Half;
|
|
|
|
Radix = FLOOR (FLOAT ("0.01") + U2 / U1);
|
|
if (Radix == E0)
|
|
printf ("Radix confirmed.\n");
|
|
else
|
|
printf ("MYSTERY: recalculated Radix = %s .\n", Radix.str());
|
|
TstCond (Defect, Radix <= Eight + Eight,
|
|
"Radix is too big: roundoff problems");
|
|
TstCond (Flaw, (Radix == Two) || (Radix == 10)
|
|
|| (Radix == One), "Radix is not as good as 2 or 10");
|
|
/*=============================================*/
|
|
Milestone = 20;
|
|
/*=============================================*/
|
|
TstCond (Failure, F9 - Half < Half,
|
|
"(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
|
|
X = F9;
|
|
I = 1;
|
|
Y = X - Half;
|
|
Z = Y - Half;
|
|
TstCond (Failure, (X != One)
|
|
|| (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
|
|
X = One + U2;
|
|
I = 0;
|
|
/*=============================================*/
|
|
Milestone = 25;
|
|
/*=============================================*/
|
|
/*... BMinusU2 = nextafter(Radix, 0) */
|
|
BMinusU2 = Radix - One;
|
|
BMinusU2 = (BMinusU2 - U2) + One;
|
|
/* Purify Integers */
|
|
if (Radix != One)
|
|
{
|
|
X = -TwoForty * LOG (U1) / LOG (Radix);
|
|
Y = FLOOR (Half + X);
|
|
if (FABS (X - Y) * Four < One)
|
|
X = Y;
|
|
Precision = X / TwoForty;
|
|
Y = FLOOR (Half + Precision);
|
|
if (FABS (Precision - Y) * TwoForty < Half)
|
|
Precision = Y;
|
|
}
|
|
if ((Precision != FLOOR (Precision)) || (Radix == One))
|
|
{
|
|
printf ("Precision cannot be characterized by an Integer number\n");
|
|
printf
|
|
("of significant digits but, by itself, this is a minor flaw.\n");
|
|
}
|
|
if (Radix == One)
|
|
printf
|
|
("logarithmic encoding has precision characterized solely by U1.\n");
|
|
else
|
|
printf ("The number of significant digits of the Radix is %s .\n",
|
|
Precision.str());
|
|
TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
|
|
"Precision worse than 5 decimal figures ");
|
|
/*=============================================*/
|
|
Milestone = 30;
|
|
/*=============================================*/
|
|
/* Test for extra-precise subexpressions */
|
|
X = FABS (((Four / Three - One) - One / Four) * Three - One / Four);
|
|
do
|
|
{
|
|
Z2 = X;
|
|
X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
|
|
}
|
|
while (!((Z2 <= X) || (X <= Zero)));
|
|
X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four);
|
|
do
|
|
{
|
|
Z1 = Z;
|
|
Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
|
|
+ One / Two)) + One / Two;
|
|
}
|
|
while (!((Z1 <= Z) || (Z <= Zero)));
|
|
do
|
|
{
|
|
do
|
|
{
|
|
Y1 = Y;
|
|
Y =
|
|
(Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half)) +
|
|
Half;
|
|
}
|
|
while (!((Y1 <= Y) || (Y <= Zero)));
|
|
X1 = X;
|
|
X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
|
|
}
|
|
while (!((X1 <= X) || (X <= Zero)));
|
|
if ((X1 != Y1) || (X1 != Z1))
|
|
{
|
|
BadCond (Serious, "Disagreements among the values X1, Y1, Z1,\n");
|
|
printf ("respectively %s, %s, %s,\n", X1.str(), Y1.str(), Z1.str());
|
|
printf ("are symptoms of inconsistencies introduced\n");
|
|
printf ("by extra-precise evaluation of arithmetic subexpressions.\n");
|
|
notify ("Possibly some part of this");
|
|
if ((X1 == U1) || (Y1 == U1) || (Z1 == U1))
|
|
printf ("That feature is not tested further by this program.\n");
|
|
}
|
|
else
|
|
{
|
|
if ((Z1 != U1) || (Z2 != U2))
|
|
{
|
|
if ((Z1 >= U1) || (Z2 >= U2))
|
|
{
|
|
BadCond (Failure, "");
|
|
notify ("Precision");
|
|
printf ("\tU1 = %s, Z1 - U1 = %s\n", U1.str(), (Z1 - U1).str());
|
|
printf ("\tU2 = %s, Z2 - U2 = %s\n", U2.str(), (Z2 - U2).str());
|
|
}
|
|
else
|
|
{
|
|
if ((Z1 <= Zero) || (Z2 <= Zero))
|
|
{
|
|
printf ("Because of unusual Radix = %s", Radix.str());
|
|
printf (", or exact rational arithmetic a result\n");
|
|
printf ("Z1 = %s, or Z2 = %s ", Z1.str(), Z2.str());
|
|
notify ("of an\nextra-precision");
|
|
}
|
|
if (Z1 != Z2 || Z1 > Zero)
|
|
{
|
|
X = Z1 / U1;
|
|
Y = Z2 / U2;
|
|
if (Y > X)
|
|
X = Y;
|
|
Q = -LOG (X);
|
|
printf ("Some subexpressions appear to be calculated "
|
|
"extra precisely\n");
|
|
printf ("with about %s extra B-digits, i.e.\n",
|
|
(Q / LOG (Radix)).str());
|
|
printf ("roughly %s extra significant decimals.\n",
|
|
(Q / LOG (FLOAT (10))).str());
|
|
}
|
|
printf
|
|
("That feature is not tested further by this program.\n");
|
|
}
|
|
}
|
|
}
|
|
Pause ();
|
|
/*=============================================*/
|
|
Milestone = 35;
|
|
/*=============================================*/
|
|
if (Radix >= Two)
|
|
{
|
|
X = W / (Radix * Radix);
|
|
Y = X + One;
|
|
Z = Y - X;
|
|
T = Z + U2;
|
|
X = T - Z;
|
|
TstCond (Failure, X == U2,
|
|
"Subtraction is not normalized X=Y,X+Z != Y+Z!");
|
|
if (X == U2)
|
|
printf ("Subtraction appears to be normalized, as it should be.");
|
|
}
|
|
printf ("\nChecking for guard digit in *, /, and -.\n");
|
|
Y = F9 * One;
|
|
Z = One * F9;
|
|
X = F9 - Half;
|
|
Y = (Y - Half) - X;
|
|
Z = (Z - Half) - X;
|
|
X = One + U2;
|
|
T = X * Radix;
|
|
R = Radix * X;
|
|
X = T - Radix;
|
|
X = X - Radix * U2;
|
|
T = R - Radix;
|
|
T = T - Radix * U2;
|
|
X = X * (Radix - One);
|
|
T = T * (Radix - One);
|
|
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero))
|
|
GMult = Yes;
|
|
else
|
|
{
|
|
GMult = No;
|
|
TstCond (Serious, false, "* lacks a Guard Digit, so 1*X != X");
|
|
}
|
|
Z = Radix * U2;
|
|
X = One + Z;
|
|
Y = FABS ((X + Z) - X * X) - U2;
|
|
X = One - U2;
|
|
Z = FABS ((X - U2) - X * X) - U1;
|
|
TstCond (Failure, (Y <= Zero)
|
|
&& (Z <= Zero), "* gets too many final digits wrong.\n");
|
|
Y = One - U2;
|
|
X = One + U2;
|
|
Z = One / Y;
|
|
Y = Z - X;
|
|
X = One / Three;
|
|
Z = Three / Nine;
|
|
X = X - Z;
|
|
T = Nine / TwentySeven;
|
|
Z = Z - T;
|
|
TstCond (Defect, X == Zero && Y == Zero && Z == Zero,
|
|
"Division lacks a Guard Digit, so error can exceed 1 ulp\n"
|
|
"or 1/3 and 3/9 and 9/27 may disagree");
|
|
Y = F9 / One;
|
|
X = F9 - Half;
|
|
Y = (Y - Half) - X;
|
|
X = One + U2;
|
|
T = X / One;
|
|
X = T - X;
|
|
if ((X == Zero) && (Y == Zero) && (Z == Zero))
|
|
GDiv = Yes;
|
|
else
|
|
{
|
|
GDiv = No;
|
|
TstCond (Serious, false, "Division lacks a Guard Digit, so X/1 != X");
|
|
}
|
|
X = One / (One + U2);
|
|
Y = X - Half - Half;
|
|
TstCond (Serious, Y < Zero, "Computed value of 1/1.000..1 >= 1");
|
|
X = One - U2;
|
|
Y = One + Radix * U2;
|
|
Z = X * Radix;
|
|
T = Y * Radix;
|
|
R = Z / Radix;
|
|
StickyBit = T / Radix;
|
|
X = R - X;
|
|
Y = StickyBit - Y;
|
|
TstCond (Failure, X == Zero && Y == Zero,
|
|
"* and/or / gets too many last digits wrong");
|
|
Y = One - U1;
|
|
X = One - F9;
|
|
Y = One - Y;
|
|
T = Radix - U2;
|
|
Z = Radix - BMinusU2;
|
|
T = Radix - T;
|
|
if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2))
|
|
GAddSub = Yes;
|
|
else
|
|
{
|
|
GAddSub = No;
|
|
TstCond (Serious, false,
|
|
"- lacks Guard Digit, so cancellation is obscured");
|
|
}
|
|
if (F9 != One && F9 - One >= Zero)
|
|
{
|
|
BadCond (Serious, "comparison alleges (1-U1) < 1 although\n");
|
|
printf (" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
|
|
printf (" such precautions against division by zero as\n");
|
|
printf (" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
|
|
}
|
|
if (GMult == Yes && GDiv == Yes && GAddSub == Yes)
|
|
printf
|
|
(" *, /, and - appear to have guard digits, as they should.\n");
|
|
/*=============================================*/
|
|
Milestone = 40;
|
|
/*=============================================*/
|
|
Pause ();
|
|
printf ("Checking rounding on multiply, divide and add/subtract.\n");
|
|
RMult = Other;
|
|
RDiv = Other;
|
|
RAddSub = Other;
|
|
RadixD2 = Radix / Two;
|
|
A1 = Two;
|
|
Done = false;
|
|
do
|
|
{
|
|
AInvrse = Radix;
|
|
do
|
|
{
|
|
X = AInvrse;
|
|
AInvrse = AInvrse / A1;
|
|
}
|
|
while (!(FLOOR (AInvrse) != AInvrse));
|
|
Done = (X == One) || (A1 > Three);
|
|
if (!Done)
|
|
A1 = Nine + One;
|
|
}
|
|
while (!(Done));
|
|
if (X == One)
|
|
A1 = Radix;
|
|
AInvrse = One / A1;
|
|
X = A1;
|
|
Y = AInvrse;
|
|
Done = false;
|
|
do
|
|
{
|
|
Z = X * Y - Half;
|
|
TstCond (Failure, Z == Half, "X * (1/X) differs from 1");
|
|
Done = X == Radix;
|
|
X = Radix;
|
|
Y = One / X;
|
|
}
|
|
while (!(Done));
|
|
Y2 = One + U2;
|
|
Y1 = One - U2;
|
|
X = OneAndHalf - U2;
|
|
Y = OneAndHalf + U2;
|
|
Z = (X - U2) * Y2;
|
|
T = Y * Y1;
|
|
Z = Z - X;
|
|
T = T - X;
|
|
X = X * Y2;
|
|
Y = (Y + U2) * Y1;
|
|
X = X - OneAndHalf;
|
|
Y = Y - OneAndHalf;
|
|
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero))
|
|
{
|
|
X = (OneAndHalf + U2) * Y2;
|
|
Y = OneAndHalf - U2 - U2;
|
|
Z = OneAndHalf + U2 + U2;
|
|
T = (OneAndHalf - U2) * Y1;
|
|
X = X - (Z + U2);
|
|
StickyBit = Y * Y1;
|
|
S = Z * Y2;
|
|
T = T - Y;
|
|
Y = (U2 - Y) + StickyBit;
|
|
Z = S - (Z + U2 + U2);
|
|
StickyBit = (Y2 + U2) * Y1;
|
|
Y1 = Y2 * Y1;
|
|
StickyBit = StickyBit - Y2;
|
|
Y1 = Y1 - Half;
|
|
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
|
|
&& (StickyBit == Zero) && (Y1 == Half))
|
|
{
|
|
RMult = Rounded;
|
|
printf ("Multiplication appears to round correctly.\n");
|
|
}
|
|
else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
|
|
&& (T < Zero) && (StickyBit + U2 == Zero) && (Y1 < Half))
|
|
{
|
|
RMult = Chopped;
|
|
printf ("Multiplication appears to chop.\n");
|
|
}
|
|
else
|
|
printf ("* is neither chopped nor correctly rounded.\n");
|
|
if ((RMult == Rounded) && (GMult == No))
|
|
notify ("Multiplication");
|
|
}
|
|
else
|
|
printf ("* is neither chopped nor correctly rounded.\n");
|
|
/*=============================================*/
|
|
Milestone = 45;
|
|
/*=============================================*/
|
|
Y2 = One + U2;
|
|
Y1 = One - U2;
|
|
Z = OneAndHalf + U2 + U2;
|
|
X = Z / Y2;
|
|
T = OneAndHalf - U2 - U2;
|
|
Y = (T - U2) / Y1;
|
|
Z = (Z + U2) / Y2;
|
|
X = X - OneAndHalf;
|
|
Y = Y - T;
|
|
T = T / Y1;
|
|
Z = Z - (OneAndHalf + U2);
|
|
T = (U2 - OneAndHalf) + T;
|
|
if (!((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero)))
|
|
{
|
|
X = OneAndHalf / Y2;
|
|
Y = OneAndHalf - U2;
|
|
Z = OneAndHalf + U2;
|
|
X = X - Y;
|
|
T = OneAndHalf / Y1;
|
|
Y = Y / Y1;
|
|
T = T - (Z + U2);
|
|
Y = Y - Z;
|
|
Z = Z / Y2;
|
|
Y1 = (Y2 + U2) / Y2;
|
|
Z = Z - OneAndHalf;
|
|
Y2 = Y1 - Y2;
|
|
Y1 = (F9 - U1) / F9;
|
|
if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
|
|
&& (Y2 == Zero) && (Y2 == Zero) && (Y1 - Half == F9 - Half))
|
|
{
|
|
RDiv = Rounded;
|
|
printf ("Division appears to round correctly.\n");
|
|
if (GDiv == No)
|
|
notify ("Division");
|
|
}
|
|
else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
|
|
&& (Y2 < Zero) && (Y1 - Half < F9 - Half))
|
|
{
|
|
RDiv = Chopped;
|
|
printf ("Division appears to chop.\n");
|
|
}
|
|
}
|
|
if (RDiv == Other)
|
|
printf ("/ is neither chopped nor correctly rounded.\n");
|
|
BInvrse = One / Radix;
|
|
TstCond (Failure, (BInvrse * Radix - Half == Half),
|
|
"Radix * ( 1 / Radix ) differs from 1");
|
|
/*=============================================*/
|
|
Milestone = 50;
|
|
/*=============================================*/
|
|
TstCond (Failure, ((F9 + U1) - Half == Half)
|
|
&& ((BMinusU2 + U2) - One == Radix - One),
|
|
"Incomplete carry-propagation in Addition");
|
|
X = One - U1 * U1;
|
|
Y = One + U2 * (One - U2);
|
|
Z = F9 - Half;
|
|
X = (X - Half) - Z;
|
|
Y = Y - One;
|
|
if ((X == Zero) && (Y == Zero))
|
|
{
|
|
RAddSub = Chopped;
|
|
printf ("Add/Subtract appears to be chopped.\n");
|
|
}
|
|
if (GAddSub == Yes)
|
|
{
|
|
X = (Half + U2) * U2;
|
|
Y = (Half - U2) * U2;
|
|
X = One + X;
|
|
Y = One + Y;
|
|
X = (One + U2) - X;
|
|
Y = One - Y;
|
|
if ((X == Zero) && (Y == Zero))
|
|
{
|
|
X = (Half + U2) * U1;
|
|
Y = (Half - U2) * U1;
|
|
X = One - X;
|
|
Y = One - Y;
|
|
X = F9 - X;
|
|
Y = One - Y;
|
|
if ((X == Zero) && (Y == Zero))
|
|
{
|
|
RAddSub = Rounded;
|
|
printf ("Addition/Subtraction appears to round correctly.\n");
|
|
if (GAddSub == No)
|
|
notify ("Add/Subtract");
|
|
}
|
|
else
|
|
printf ("Addition/Subtraction neither rounds nor chops.\n");
|
|
}
|
|
else
|
|
printf ("Addition/Subtraction neither rounds nor chops.\n");
|
|
}
|
|
else
|
|
printf ("Addition/Subtraction neither rounds nor chops.\n");
|
|
S = One;
|
|
X = One + Half * (One + Half);
|
|
Y = (One + U2) * Half;
|
|
Z = X - Y;
|
|
T = Y - X;
|
|
StickyBit = Z + T;
|
|
if (StickyBit != Zero)
|
|
{
|
|
S = Zero;
|
|
BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n");
|
|
}
|
|
StickyBit = Zero;
|
|
if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
|
|
&& (RMult == Rounded) && (RDiv == Rounded)
|
|
&& (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2))
|
|
{
|
|
printf ("Checking for sticky bit.\n");
|
|
X = (Half + U1) * U2;
|
|
Y = Half * U2;
|
|
Z = One + Y;
|
|
T = One + X;
|
|
if ((Z - One <= Zero) && (T - One >= U2))
|
|
{
|
|
Z = T + Y;
|
|
Y = Z - X;
|
|
if ((Z - T >= U2) && (Y - T == Zero))
|
|
{
|
|
X = (Half + U1) * U1;
|
|
Y = Half * U1;
|
|
Z = One - Y;
|
|
T = One - X;
|
|
if ((Z - One == Zero) && (T - F9 == Zero))
|
|
{
|
|
Z = (Half - U1) * U1;
|
|
T = F9 - Z;
|
|
Q = F9 - Y;
|
|
if ((T - F9 == Zero) && (F9 - U1 - Q == Zero))
|
|
{
|
|
Z = (One + U2) * OneAndHalf;
|
|
T = (OneAndHalf + U2) - Z + U2;
|
|
X = One + Half / Radix;
|
|
Y = One + Radix * U2;
|
|
Z = X * Y;
|
|
if (T == Zero && X + Radix * U2 - Z == Zero)
|
|
{
|
|
if (Radix != Two)
|
|
{
|
|
X = Two + U2;
|
|
Y = X / Two;
|
|
if ((Y - One == Zero))
|
|
StickyBit = S;
|
|
}
|
|
else
|
|
StickyBit = S;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (StickyBit == One)
|
|
printf ("Sticky bit apparently used correctly.\n");
|
|
else
|
|
printf ("Sticky bit used incorrectly or not at all.\n");
|
|
TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
|
|
RMult == Other || RDiv == Other || RAddSub == Other),
|
|
"lack(s) of guard digits or failure(s) to correctly round or chop\n\
|
|
(noted above) count as one flaw in the final tally below");
|
|
/*=============================================*/
|
|
Milestone = 60;
|
|
/*=============================================*/
|
|
printf ("\n");
|
|
printf ("Does Multiplication commute? ");
|
|
printf ("Testing on %d random pairs.\n", NoTrials);
|
|
Random9 = SQRT (FLOAT (3));
|
|
Random1 = Third;
|
|
I = 1;
|
|
do
|
|
{
|
|
X = Random ();
|
|
Y = Random ();
|
|
Z9 = Y * X;
|
|
Z = X * Y;
|
|
Z9 = Z - Z9;
|
|
I = I + 1;
|
|
}
|
|
while (!((I > NoTrials) || (Z9 != Zero)));
|
|
if (I == NoTrials)
|
|
{
|
|
Random1 = One + Half / Three;
|
|
Random2 = (U2 + U1) + One;
|
|
Z = Random1 * Random2;
|
|
Y = Random2 * Random1;
|
|
Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
|
|
Three) * ((U2 + U1) +
|
|
One);
|
|
}
|
|
if (!((I == NoTrials) || (Z9 == Zero)))
|
|
BadCond (Defect, "X * Y == Y * X trial fails.\n");
|
|
else
|
|
printf (" No failures found in %d integer pairs.\n", NoTrials);
|
|
/*=============================================*/
|
|
Milestone = 70;
|
|
/*=============================================*/
|
|
printf ("\nRunning test of square root(x).\n");
|
|
TstCond (Failure, (Zero == SQRT (Zero))
|
|
&& (-Zero == SQRT (-Zero))
|
|
&& (One == SQRT (One)), "Square root of 0.0, -0.0 or 1.0 wrong");
|
|
MinSqEr = Zero;
|
|
MaxSqEr = Zero;
|
|
J = Zero;
|
|
X = Radix;
|
|
OneUlp = U2;
|
|
SqXMinX (Serious);
|
|
X = BInvrse;
|
|
OneUlp = BInvrse * U1;
|
|
SqXMinX (Serious);
|
|
X = U1;
|
|
OneUlp = U1 * U1;
|
|
SqXMinX (Serious);
|
|
if (J != Zero)
|
|
Pause ();
|
|
printf ("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
|
|
J = Zero;
|
|
X = Two;
|
|
Y = Radix;
|
|
if ((Radix != One))
|
|
do
|
|
{
|
|
X = Y;
|
|
Y = Radix * Y;
|
|
}
|
|
while (!((Y - X >= NoTrials)));
|
|
OneUlp = X * U2;
|
|
I = 1;
|
|
while (I <= NoTrials)
|
|
{
|
|
X = X + One;
|
|
SqXMinX (Defect);
|
|
if (J > Zero)
|
|
break;
|
|
I = I + 1;
|
|
}
|
|
printf ("Test for sqrt monotonicity.\n");
|
|
I = -1;
|
|
X = BMinusU2;
|
|
Y = Radix;
|
|
Z = Radix + Radix * U2;
|
|
NotMonot = false;
|
|
Monot = false;
|
|
while (!(NotMonot || Monot))
|
|
{
|
|
I = I + 1;
|
|
X = SQRT (X);
|
|
Q = SQRT (Y);
|
|
Z = SQRT (Z);
|
|
if ((X > Q) || (Q > Z))
|
|
NotMonot = true;
|
|
else
|
|
{
|
|
Q = FLOOR (Q + Half);
|
|
if (!(I > 0 || Radix == Q * Q))
|
|
Monot = true;
|
|
else if (I > 0)
|
|
{
|
|
if (I > 1)
|
|
Monot = true;
|
|
else
|
|
{
|
|
Y = Y * BInvrse;
|
|
X = Y - U1;
|
|
Z = Y + U1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Y = Q;
|
|
X = Y - U2;
|
|
Z = Y + U2;
|
|
}
|
|
}
|
|
}
|
|
if (Monot)
|
|
printf ("sqrt has passed a test for Monotonicity.\n");
|
|
else
|
|
{
|
|
BadCond (Defect, "");
|
|
printf ("sqrt(X) is non-monotonic for X near %s .\n", Y.str());
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 110;
|
|
/*=============================================*/
|
|
printf ("Seeking Underflow thresholds UfThold and E0.\n");
|
|
D = U1;
|
|
if (Precision != FLOOR (Precision))
|
|
{
|
|
D = BInvrse;
|
|
X = Precision;
|
|
do
|
|
{
|
|
D = D * BInvrse;
|
|
X = X - One;
|
|
}
|
|
while (X > Zero);
|
|
}
|
|
Y = One;
|
|
Z = D;
|
|
/* ... D is power of 1/Radix < 1. */
|
|
do
|
|
{
|
|
C = Y;
|
|
Y = Z;
|
|
Z = Y * Y;
|
|
}
|
|
while ((Y > Z) && (Z + Z > Z));
|
|
Y = C;
|
|
Z = Y * D;
|
|
do
|
|
{
|
|
C = Y;
|
|
Y = Z;
|
|
Z = Y * D;
|
|
}
|
|
while ((Y > Z) && (Z + Z > Z));
|
|
if (Radix < Two)
|
|
HInvrse = Two;
|
|
else
|
|
HInvrse = Radix;
|
|
H = One / HInvrse;
|
|
/* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
|
|
CInvrse = One / C;
|
|
E0 = C;
|
|
Z = E0 * H;
|
|
/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
|
|
do
|
|
{
|
|
Y = E0;
|
|
E0 = Z;
|
|
Z = E0 * H;
|
|
}
|
|
while ((E0 > Z) && (Z + Z > Z));
|
|
UfThold = E0;
|
|
E1 = Zero;
|
|
Q = Zero;
|
|
E9 = U2;
|
|
S = One + E9;
|
|
D = C * S;
|
|
if (D <= C)
|
|
{
|
|
E9 = Radix * U2;
|
|
S = One + E9;
|
|
D = C * S;
|
|
if (D <= C)
|
|
{
|
|
BadCond (Failure,
|
|
"multiplication gets too many last digits wrong.\n");
|
|
Underflow = E0;
|
|
Y1 = Zero;
|
|
PseudoZero = Z;
|
|
Pause ();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Underflow = D;
|
|
PseudoZero = Underflow * H;
|
|
UfThold = Zero;
|
|
do
|
|
{
|
|
Y1 = Underflow;
|
|
Underflow = PseudoZero;
|
|
if (E1 + E1 <= E1)
|
|
{
|
|
Y2 = Underflow * HInvrse;
|
|
E1 = FABS (Y1 - Y2);
|
|
Q = Y1;
|
|
if ((UfThold == Zero) && (Y1 != Y2))
|
|
UfThold = Y1;
|
|
}
|
|
PseudoZero = PseudoZero * H;
|
|
}
|
|
while ((Underflow > PseudoZero)
|
|
&& (PseudoZero + PseudoZero > PseudoZero));
|
|
}
|
|
/* Comment line 4530 .. 4560 */
|
|
if (PseudoZero != Zero)
|
|
{
|
|
printf ("\n");
|
|
Z = PseudoZero;
|
|
/* ... Test PseudoZero for "phoney- zero" violates */
|
|
/* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
|
|
... */
|
|
if (PseudoZero <= Zero)
|
|
{
|
|
BadCond (Failure, "Positive expressions can underflow to an\n");
|
|
printf ("allegedly negative value\n");
|
|
printf ("PseudoZero that prints out as: %s .\n", PseudoZero.str());
|
|
X = -PseudoZero;
|
|
if (X <= Zero)
|
|
{
|
|
printf ("But -PseudoZero, which should be\n");
|
|
printf ("positive, isn't; it prints out as %s .\n", X.str());
|
|
}
|
|
}
|
|
else
|
|
{
|
|
BadCond (Flaw, "Underflow can stick at an allegedly positive\n");
|
|
printf ("value PseudoZero that prints out as %s .\n",
|
|
PseudoZero.str());
|
|
}
|
|
TstPtUf ();
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 120;
|
|
/*=============================================*/
|
|
if (CInvrse * Y > CInvrse * Y1)
|
|
{
|
|
S = H * S;
|
|
E0 = Underflow;
|
|
}
|
|
if (!((E1 == Zero) || (E1 == E0)))
|
|
{
|
|
BadCond (Defect, "");
|
|
if (E1 < E0)
|
|
{
|
|
printf ("Products underflow at a higher");
|
|
printf (" threshold than differences.\n");
|
|
if (PseudoZero == Zero)
|
|
E0 = E1;
|
|
}
|
|
else
|
|
{
|
|
printf ("Difference underflows at a higher");
|
|
printf (" threshold than products.\n");
|
|
}
|
|
}
|
|
printf ("Smallest strictly positive number found is E0 = %s .\n", E0.str());
|
|
Z = E0;
|
|
TstPtUf ();
|
|
Underflow = E0;
|
|
if (N == 1)
|
|
Underflow = Y;
|
|
I = 4;
|
|
if (E1 == Zero)
|
|
I = 3;
|
|
if (UfThold == Zero)
|
|
I = I - 2;
|
|
UfNGrad = true;
|
|
switch (I)
|
|
{
|
|
case 1:
|
|
UfThold = Underflow;
|
|
if ((CInvrse * Q) != ((CInvrse * Y) * S))
|
|
{
|
|
UfThold = Y;
|
|
BadCond (Failure, "Either accuracy deteriorates as numbers\n");
|
|
printf ("approach a threshold = %s\n", UfThold.str());
|
|
printf (" coming down from %s\n", C.str());
|
|
printf
|
|
(" or else multiplication gets too many last digits wrong.\n");
|
|
}
|
|
Pause ();
|
|
break;
|
|
|
|
case 2:
|
|
BadCond (Failure,
|
|
"Underflow confuses Comparison, which alleges that\n");
|
|
printf ("Q == Y while denying that |Q - Y| == 0; these values\n");
|
|
printf ("print out as Q = %s, Y = %s .\n", Q.str(), Y2.str());
|
|
printf ("|Q - Y| = %s .\n", FABS (Q - Y2).str());
|
|
UfThold = Q;
|
|
break;
|
|
|
|
case 3:
|
|
X = X;
|
|
break;
|
|
|
|
case 4:
|
|
if ((Q == UfThold) && (E1 == E0) && (FABS (UfThold - E1 / E9) <= E1))
|
|
{
|
|
UfNGrad = false;
|
|
printf ("Underflow is gradual; it incurs Absolute Error =\n");
|
|
printf ("(roundoff in UfThold) < E0.\n");
|
|
Y = E0 * CInvrse;
|
|
Y = Y * (OneAndHalf + U2);
|
|
X = CInvrse * (One + U2);
|
|
Y = Y / X;
|
|
IEEE = (Y == E0);
|
|
}
|
|
}
|
|
if (UfNGrad)
|
|
{
|
|
printf ("\n");
|
|
if (setjmp (ovfl_buf))
|
|
{
|
|
printf ("Underflow / UfThold failed!\n");
|
|
R = H + H;
|
|
}
|
|
else
|
|
R = SQRT (Underflow / UfThold);
|
|
if (R <= H)
|
|
{
|
|
Z = R * UfThold;
|
|
X = Z * (One + R * H * (One + H));
|
|
}
|
|
else
|
|
{
|
|
Z = UfThold;
|
|
X = Z * (One + H * H * (One + H));
|
|
}
|
|
if (!((X == Z) || (X - Z != Zero)))
|
|
{
|
|
BadCond (Flaw, "");
|
|
printf ("X = %s\n\tis not equal to Z = %s .\n", X.str(), Z.str());
|
|
Z9 = X - Z;
|
|
printf ("yet X - Z yields %s .\n", Z9.str());
|
|
printf (" Should this NOT signal Underflow, ");
|
|
printf ("this is a SERIOUS DEFECT\nthat causes ");
|
|
printf ("confusion when innocent statements like\n");;
|
|
printf (" if (X == Z) ... else");
|
|
printf (" ... (f(X) - f(Z)) / (X - Z) ...\n");
|
|
printf ("encounter Division by Zero although actually\n");
|
|
if (setjmp (ovfl_buf))
|
|
printf ("X / Z fails!\n");
|
|
else
|
|
printf ("X / Z = 1 + %s .\n", ((X / Z - Half) - Half).str());
|
|
}
|
|
}
|
|
printf ("The Underflow threshold is %s, below which\n", UfThold.str());
|
|
printf ("calculation may suffer larger Relative error than ");
|
|
printf ("merely roundoff.\n");
|
|
Y2 = U1 * U1;
|
|
Y = Y2 * Y2;
|
|
Y2 = Y * U1;
|
|
if (Y2 <= UfThold)
|
|
{
|
|
if (Y > E0)
|
|
{
|
|
BadCond (Defect, "");
|
|
I = 5;
|
|
}
|
|
else
|
|
{
|
|
BadCond (Serious, "");
|
|
I = 4;
|
|
}
|
|
printf ("Range is too narrow; U1^%d Underflows.\n", I);
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 130;
|
|
/*=============================================*/
|
|
Y = -FLOOR (Half - TwoForty * LOG (UfThold) / LOG (HInvrse)) / TwoForty;
|
|
Y2 = Y + Y;
|
|
printf ("Since underflow occurs below the threshold\n");
|
|
printf ("UfThold = (%s) ^ (%s)\nonly underflow ", HInvrse.str(), Y.str());
|
|
printf ("should afflict the expression\n\t(%s) ^ (%s);\n",
|
|
HInvrse.str(), Y2.str());
|
|
printf ("actually calculating yields:");
|
|
if (setjmp (ovfl_buf))
|
|
{
|
|
BadCond (Serious, "trap on underflow.\n");
|
|
}
|
|
else
|
|
{
|
|
V9 = POW (HInvrse, Y2);
|
|
printf (" %s .\n", V9.str());
|
|
if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold)))
|
|
{
|
|
BadCond (Serious, "this is not between 0 and underflow\n");
|
|
printf (" threshold = %s .\n", UfThold.str());
|
|
}
|
|
else if (!(V9 > UfThold * (One + E9)))
|
|
printf ("This computed value is O.K.\n");
|
|
else
|
|
{
|
|
BadCond (Defect, "this is not between 0 and underflow\n");
|
|
printf (" threshold = %s .\n", UfThold.str());
|
|
}
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 160;
|
|
/*=============================================*/
|
|
Pause ();
|
|
printf ("Searching for Overflow threshold:\n");
|
|
printf ("This may generate an error.\n");
|
|
Y = -CInvrse;
|
|
V9 = HInvrse * Y;
|
|
if (setjmp (ovfl_buf))
|
|
{
|
|
I = 0;
|
|
V9 = Y;
|
|
goto overflow;
|
|
}
|
|
do
|
|
{
|
|
V = Y;
|
|
Y = V9;
|
|
V9 = HInvrse * Y;
|
|
}
|
|
while (V9 < Y);
|
|
I = 1;
|
|
overflow:
|
|
Z = V9;
|
|
printf ("Can `Z = -Y' overflow?\n");
|
|
printf ("Trying it on Y = %s .\n", Y.str());
|
|
V9 = -Y;
|
|
V0 = V9;
|
|
if (V - Y == V + V0)
|
|
printf ("Seems O.K.\n");
|
|
else
|
|
{
|
|
printf ("finds a ");
|
|
BadCond (Flaw, "-(-Y) differs from Y.\n");
|
|
}
|
|
if (Z != Y)
|
|
{
|
|
BadCond (Serious, "");
|
|
printf ("overflow past %s\n\tshrinks to %s .\n", Y.str(), Z.str());
|
|
}
|
|
if (I)
|
|
{
|
|
Y = V * (HInvrse * U2 - HInvrse);
|
|
Z = Y + ((One - HInvrse) * U2) * V;
|
|
if (Z < V0)
|
|
Y = Z;
|
|
if (Y < V0)
|
|
V = Y;
|
|
if (V0 - V < V0)
|
|
V = V0;
|
|
}
|
|
else
|
|
{
|
|
V = Y * (HInvrse * U2 - HInvrse);
|
|
V = V + ((One - HInvrse) * U2) * Y;
|
|
}
|
|
printf ("Overflow threshold is V = %s .\n", V.str());
|
|
if (I)
|
|
printf ("Overflow saturates at V0 = %s .\n", V0.str());
|
|
else
|
|
printf ("There is no saturation value because "
|
|
"the system traps on overflow.\n");
|
|
V9 = V * One;
|
|
printf ("No Overflow should be signaled for V * 1 = %s\n", V9.str());
|
|
V9 = V / One;
|
|
printf (" nor for V / 1 = %s.\n", V9.str());
|
|
printf ("Any overflow signal separating this * from the one\n");
|
|
printf ("above is a DEFECT.\n");
|
|
/*=============================================*/
|
|
Milestone = 170;
|
|
/*=============================================*/
|
|
if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V))
|
|
{
|
|
BadCond (Failure, "Comparisons involving ");
|
|
printf ("+-%s, +-%s\nand +-%s are confused by Overflow.",
|
|
V.str(), V0.str(), UfThold.str());
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 175;
|
|
/*=============================================*/
|
|
printf ("\n");
|
|
for (Indx = 1; Indx <= 3; ++Indx)
|
|
{
|
|
switch (Indx)
|
|
{
|
|
case 1:
|
|
Z = UfThold;
|
|
break;
|
|
case 2:
|
|
Z = E0;
|
|
break;
|
|
case 3:
|
|
Z = PseudoZero;
|
|
break;
|
|
}
|
|
if (Z != Zero)
|
|
{
|
|
V9 = SQRT (Z);
|
|
Y = V9 * V9;
|
|
if (Y / (One - Radix * E9) < Z || Y > (One + Radix * E9) * Z)
|
|
{ /* dgh: + E9 --> * E9 */
|
|
if (V9 > U1)
|
|
BadCond (Serious, "");
|
|
else
|
|
BadCond (Defect, "");
|
|
printf ("Comparison alleges that what prints as Z = %s\n",
|
|
Z.str());
|
|
printf (" is too far from sqrt(Z) ^ 2 = %s .\n", Y.str());
|
|
}
|
|
}
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 180;
|
|
/*=============================================*/
|
|
for (Indx = 1; Indx <= 2; ++Indx)
|
|
{
|
|
if (Indx == 1)
|
|
Z = V;
|
|
else
|
|
Z = V0;
|
|
V9 = SQRT (Z);
|
|
X = (One - Radix * E9) * V9;
|
|
V9 = V9 * X;
|
|
if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z)))
|
|
{
|
|
Y = V9;
|
|
if (X < W)
|
|
BadCond (Serious, "");
|
|
else
|
|
BadCond (Defect, "");
|
|
printf ("Comparison alleges that Z = %s\n", Z.str());
|
|
printf (" is too far from sqrt(Z) ^ 2 (%s) .\n", Y.str());
|
|
}
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 190;
|
|
/*=============================================*/
|
|
Pause ();
|
|
X = UfThold * V;
|
|
Y = Radix * Radix;
|
|
if (X * Y < One || X > Y)
|
|
{
|
|
if (X * Y < U1 || X > Y / U1)
|
|
BadCond (Defect, "Badly");
|
|
else
|
|
BadCond (Flaw, "");
|
|
|
|
printf (" unbalanced range; UfThold * V = %s\n\t%s\n",
|
|
X.str(), "is too far from 1.\n");
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 200;
|
|
/*=============================================*/
|
|
for (Indx = 1; Indx <= 5; ++Indx)
|
|
{
|
|
X = F9;
|
|
switch (Indx)
|
|
{
|
|
case 2:
|
|
X = One + U2;
|
|
break;
|
|
case 3:
|
|
X = V;
|
|
break;
|
|
case 4:
|
|
X = UfThold;
|
|
break;
|
|
case 5:
|
|
X = Radix;
|
|
}
|
|
Y = X;
|
|
if (setjmp (ovfl_buf))
|
|
printf (" X / X traps when X = %s\n", X.str());
|
|
else
|
|
{
|
|
V9 = (Y / X - Half) - Half;
|
|
if (V9 == Zero)
|
|
continue;
|
|
if (V9 == -U1 && Indx < 5)
|
|
BadCond (Flaw, "");
|
|
else
|
|
BadCond (Serious, "");
|
|
printf (" X / X differs from 1 when X = %s\n", X.str());
|
|
printf (" instead, X / X - 1/2 - 1/2 = %s .\n", V9.str());
|
|
}
|
|
}
|
|
/*=============================================*/
|
|
Milestone = 210;
|
|
/*=============================================*/
|
|
MyZero = Zero;
|
|
printf ("\n");
|
|
printf ("What message and/or values does Division by Zero produce?\n");
|
|
printf (" Trying to compute 1 / 0 produces ...");
|
|
if (!setjmp (ovfl_buf))
|
|
printf (" %s .\n", (One / MyZero).str());
|
|
printf ("\n Trying to compute 0 / 0 produces ...");
|
|
if (!setjmp (ovfl_buf))
|
|
printf (" %s .\n", (Zero / MyZero).str());
|
|
/*=============================================*/
|
|
Milestone = 220;
|
|
/*=============================================*/
|
|
Pause ();
|
|
printf ("\n");
|
|
{
|
|
static const char *msg[] = {
|
|
"FAILUREs encountered =",
|
|
"SERIOUS DEFECTs discovered =",
|
|
"DEFECTs discovered =",
|
|
"FLAWs discovered ="
|
|
};
|
|
int i;
|
|
for (i = 0; i < 4; i++)
|
|
if (ErrCnt[i])
|
|
printf ("The number of %-29s %d.\n", msg[i], ErrCnt[i]);
|
|
}
|
|
printf ("\n");
|
|
if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] + ErrCnt[Flaw]) > 0)
|
|
{
|
|
if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] == 0)
|
|
&& (ErrCnt[Flaw] > 0))
|
|
{
|
|
printf ("The arithmetic diagnosed seems ");
|
|
printf ("Satisfactory though flawed.\n");
|
|
}
|
|
if ((ErrCnt[Failure] + ErrCnt[Serious] == 0) && (ErrCnt[Defect] > 0))
|
|
{
|
|
printf ("The arithmetic diagnosed may be Acceptable\n");
|
|
printf ("despite inconvenient Defects.\n");
|
|
}
|
|
if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0)
|
|
{
|
|
printf ("The arithmetic diagnosed has ");
|
|
printf ("unacceptable Serious Defects.\n");
|
|
}
|
|
if (ErrCnt[Failure] > 0)
|
|
{
|
|
printf ("Potentially fatal FAILURE may have spoiled this");
|
|
printf (" program's subsequent diagnoses.\n");
|
|
}
|
|
}
|
|
else
|
|
{
|
|
printf ("No failures, defects nor flaws have been discovered.\n");
|
|
if (!((RMult == Rounded) && (RDiv == Rounded)
|
|
&& (RAddSub == Rounded) && (RSqrt == Rounded)))
|
|
printf ("The arithmetic diagnosed seems Satisfactory.\n");
|
|
else
|
|
{
|
|
if (StickyBit >= One &&
|
|
(Radix - Two) * (Radix - Nine - One) == Zero)
|
|
{
|
|
printf ("Rounding appears to conform to ");
|
|
printf ("the proposed IEEE standard P");
|
|
if ((Radix == Two) &&
|
|
((Precision - Four * Three * Two) *
|
|
(Precision - TwentySeven - TwentySeven + One) == Zero))
|
|
printf ("754");
|
|
else
|
|
printf ("854");
|
|
if (IEEE)
|
|
printf (".\n");
|
|
else
|
|
{
|
|
printf (",\nexcept for possibly Double Rounding");
|
|
printf (" during Gradual Underflow.\n");
|
|
}
|
|
}
|
|
printf ("The arithmetic diagnosed appears to be Excellent!\n");
|
|
}
|
|
}
|
|
printf ("END OF TEST.\n");
|
|
return 0;
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
Paranoia<FLOAT>::Sign (FLOAT X)
|
|
{
|
|
return X >= FLOAT (long (0)) ? 1 : -1;
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::Pause ()
|
|
{
|
|
if (do_pause)
|
|
{
|
|
fputs ("Press return...", stdout);
|
|
fflush (stdout);
|
|
getchar();
|
|
}
|
|
printf ("\nDiagnosis resumes after milestone Number %d", Milestone);
|
|
printf (" Page: %d\n\n", PageNo);
|
|
++Milestone;
|
|
++PageNo;
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::TstCond (int K, int Valid, const char *T)
|
|
{
|
|
if (!Valid)
|
|
{
|
|
BadCond (K, T);
|
|
printf (".\n");
|
|
}
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::BadCond (int K, const char *T)
|
|
{
|
|
static const char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
|
|
|
|
ErrCnt[K] = ErrCnt[K] + 1;
|
|
printf ("%s: %s", msg[K], T);
|
|
}
|
|
|
|
/* Random computes
|
|
X = (Random1 + Random9)^5
|
|
Random1 = X - FLOOR(X) + 0.000005 * X;
|
|
and returns the new value of Random1. */
|
|
|
|
template<typename FLOAT>
|
|
FLOAT
|
|
Paranoia<FLOAT>::Random ()
|
|
{
|
|
FLOAT X, Y;
|
|
|
|
X = Random1 + Random9;
|
|
Y = X * X;
|
|
Y = Y * Y;
|
|
X = X * Y;
|
|
Y = X - FLOOR (X);
|
|
Random1 = Y + X * FLOAT ("0.000005");
|
|
return (Random1);
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::SqXMinX (int ErrKind)
|
|
{
|
|
FLOAT XA, XB;
|
|
|
|
XB = X * BInvrse;
|
|
XA = X - XB;
|
|
SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp;
|
|
if (SqEr != Zero)
|
|
{
|
|
if (SqEr < MinSqEr)
|
|
MinSqEr = SqEr;
|
|
if (SqEr > MaxSqEr)
|
|
MaxSqEr = SqEr;
|
|
J = J + 1;
|
|
BadCond (ErrKind, "\n");
|
|
printf ("sqrt(%s) - %s = %s\n", (X * X).str(), X.str(),
|
|
(OneUlp * SqEr).str());
|
|
printf ("\tinstead of correct value 0 .\n");
|
|
}
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::NewD ()
|
|
{
|
|
X = Z1 * Q;
|
|
X = FLOOR (Half - X / Radix) * Radix + X;
|
|
Q = (Q - X * Z) / Radix + X * X * (D / Radix);
|
|
Z = Z - Two * X * D;
|
|
if (Z <= Zero)
|
|
{
|
|
Z = -Z;
|
|
Z1 = -Z1;
|
|
}
|
|
D = Radix * D;
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::SR3750 ()
|
|
{
|
|
if (!((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2)))
|
|
{
|
|
I = I + 1;
|
|
X2 = SQRT (X * D);
|
|
Y2 = (X2 - Z2) - (Y - Z2);
|
|
X2 = X8 / (Y - Half);
|
|
X2 = X2 - Half * X2 * X2;
|
|
SqEr = (Y2 + Half) + (Half - X2);
|
|
if (SqEr < MinSqEr)
|
|
MinSqEr = SqEr;
|
|
SqEr = Y2 - X2;
|
|
if (SqEr > MaxSqEr)
|
|
MaxSqEr = SqEr;
|
|
}
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::IsYeqX ()
|
|
{
|
|
if (Y != X)
|
|
{
|
|
if (N <= 0)
|
|
{
|
|
if (Z == Zero && Q <= Zero)
|
|
printf ("WARNING: computing\n");
|
|
else
|
|
BadCond (Defect, "computing\n");
|
|
printf ("\t(%s) ^ (%s)\n", Z.str(), Q.str());
|
|
printf ("\tyielded %s;\n", Y.str());
|
|
printf ("\twhich compared unequal to correct %s ;\n", X.str());
|
|
printf ("\t\tthey differ by %s .\n", (Y - X).str());
|
|
}
|
|
N = N + 1; /* ... count discrepancies. */
|
|
}
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::PrintIfNPositive ()
|
|
{
|
|
if (N > 0)
|
|
printf ("Similar discrepancies have occurred %d times.\n", N);
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::TstPtUf ()
|
|
{
|
|
N = 0;
|
|
if (Z != Zero)
|
|
{
|
|
printf ("Since comparison denies Z = 0, evaluating ");
|
|
printf ("(Z + Z) / Z should be safe.\n");
|
|
if (setjmp (ovfl_buf))
|
|
goto very_serious;
|
|
Q9 = (Z + Z) / Z;
|
|
printf ("What the machine gets for (Z + Z) / Z is %s .\n", Q9.str());
|
|
if (FABS (Q9 - Two) < Radix * U2)
|
|
{
|
|
printf ("This is O.K., provided Over/Underflow");
|
|
printf (" has NOT just been signaled.\n");
|
|
}
|
|
else
|
|
{
|
|
if ((Q9 < One) || (Q9 > Two))
|
|
{
|
|
very_serious:
|
|
N = 1;
|
|
ErrCnt[Serious] = ErrCnt[Serious] + 1;
|
|
printf ("This is a VERY SERIOUS DEFECT!\n");
|
|
}
|
|
else
|
|
{
|
|
N = 1;
|
|
ErrCnt[Defect] = ErrCnt[Defect] + 1;
|
|
printf ("This is a DEFECT!\n");
|
|
}
|
|
}
|
|
V9 = Z * One;
|
|
Random1 = V9;
|
|
V9 = One * Z;
|
|
Random2 = V9;
|
|
V9 = Z / One;
|
|
if ((Z == Random1) && (Z == Random2) && (Z == V9))
|
|
{
|
|
if (N > 0)
|
|
Pause ();
|
|
}
|
|
else
|
|
{
|
|
N = 1;
|
|
BadCond (Defect, "What prints as Z = ");
|
|
printf ("%s\n\tcompares different from ", Z.str());
|
|
if (Z != Random1)
|
|
printf ("Z * 1 = %s ", Random1.str());
|
|
if (!((Z == Random2) || (Random2 == Random1)))
|
|
printf ("1 * Z == %s\n", Random2.str());
|
|
if (!(Z == V9))
|
|
printf ("Z / 1 = %s\n", V9.str());
|
|
if (Random2 != Random1)
|
|
{
|
|
ErrCnt[Defect] = ErrCnt[Defect] + 1;
|
|
BadCond (Defect, "Multiplication does not commute!\n");
|
|
printf ("\tComparison alleges that 1 * Z = %s\n", Random2.str());
|
|
printf ("\tdiffers from Z * 1 = %s\n", Random1.str());
|
|
}
|
|
Pause ();
|
|
}
|
|
}
|
|
}
|
|
|
|
template<typename FLOAT>
|
|
void
|
|
Paranoia<FLOAT>::notify (const char *s)
|
|
{
|
|
printf ("%s test appears to be inconsistent...\n", s);
|
|
printf (" PLEASE NOTIFY KARPINKSI!\n");
|
|
}
|
|
|
|
/* ====================================================================== */
|
|
|
|
int main(int ac, char **av)
|
|
{
|
|
setbuf(stdout, NULL);
|
|
setbuf(stderr, NULL);
|
|
|
|
while (1)
|
|
switch (getopt (ac, av, "pvg:fdl"))
|
|
{
|
|
case -1:
|
|
return 0;
|
|
case 'p':
|
|
do_pause = true;
|
|
break;
|
|
case 'v':
|
|
verbose = true;
|
|
break;
|
|
case 'g':
|
|
{
|
|
static const struct {
|
|
const char *name;
|
|
const struct real_format *fmt;
|
|
} fmts[] = {
|
|
#define F(x) { #x, &x##_format }
|
|
F(ieee_single),
|
|
F(ieee_double),
|
|
F(ieee_extended_motorola),
|
|
F(ieee_extended_intel_96),
|
|
F(ieee_extended_intel_128),
|
|
F(ibm_extended),
|
|
F(ieee_quad),
|
|
F(vax_f),
|
|
F(vax_d),
|
|
F(vax_g),
|
|
F(i370_single),
|
|
F(i370_double),
|
|
F(c4x_single),
|
|
F(c4x_extended),
|
|
F(real_internal),
|
|
#undef F
|
|
};
|
|
|
|
int i, n = sizeof (fmts)/sizeof(*fmts);
|
|
|
|
for (i = 0; i < n; ++i)
|
|
if (strcmp (fmts[i].name, optarg) == 0)
|
|
break;
|
|
|
|
if (i == n)
|
|
{
|
|
printf ("Unknown implementation \"%s\"; "
|
|
"available implementations:\n", optarg);
|
|
for (i = 0; i < n; ++i)
|
|
printf ("\t%s\n", fmts[i].name);
|
|
return 1;
|
|
}
|
|
|
|
// We cheat and use the same mode all the time, but vary
|
|
// the format used for that mode.
|
|
real_format_for_mode[int(real_c_float::MODE) - int(QFmode)]
|
|
= fmts[i].fmt;
|
|
|
|
Paranoia<real_c_float>().main();
|
|
break;
|
|
}
|
|
|
|
case 'f':
|
|
Paranoia < native_float<float> >().main();
|
|
break;
|
|
case 'd':
|
|
Paranoia < native_float<double> >().main();
|
|
break;
|
|
case 'l':
|
|
#ifndef NO_LONG_DOUBLE
|
|
Paranoia < native_float<long double> >().main();
|
|
#endif
|
|
break;
|
|
|
|
case '?':
|
|
puts ("-p\tpause between pages");
|
|
puts ("-g<FMT>\treal.c implementation FMT");
|
|
puts ("-f\tnative float");
|
|
puts ("-d\tnative double");
|
|
puts ("-l\tnative long double");
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
/* GCC stuff referenced by real.o. */
|
|
|
|
extern "C" void
|
|
fancy_abort ()
|
|
{
|
|
abort ();
|
|
}
|
|
|
|
int target_flags = 0;
|
|
|
|
extern "C" int
|
|
floor_log2_wide (unsigned HOST_WIDE_INT x)
|
|
{
|
|
int log = -1;
|
|
while (x != 0)
|
|
log++,
|
|
x >>= 1;
|
|
return log;
|
|
}
|