contrib/paranoia.cc

   1 /*      A C version of Kahan's Floating Point Test "Paranoia"
   2
   3 Thos Sumner, UCSF, Feb. 1985
   4 David Gay, BTL, Jan. 1986
   5
   6 This is a rewrite from the Pascal version by
   7
   8 B. A. Wichmann, 18 Jan. 1985
   9
  10 (and does NOT exhibit good C programming style).
  11
  12 Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
  13
  14 (C) Apr 19 1983 in BASIC version by:
  15 Professor W. M. Kahan,
  16 567 Evans Hall
  17 Electrical Engineering & Computer Science Dept.
  18 University of California
  19 Berkeley, California 94720
  20 USA
  21
  22 converted to Pascal by:
  23 B. A. Wichmann
  24 National Physical Laboratory
  25 Teddington Middx
  26 TW11 OLW
  27 UK
  28
  29 converted to C by:
  30
  31 David M. Gay            and     Thos Sumner
  32 AT&T Bell Labs                  Computer Center, Rm. U-76
  33 600 Mountain Avenue             University of California
  34 Murray Hill, NJ 07974           San Francisco, CA 94143
  35 USA                             USA
  36
  37 with simultaneous corrections to the Pascal source (reflected
  38 in the Pascal source available over netlib).
  39 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
  40
  41 Reports of results on various systems from all the versions
  42 of Paranoia are being collected by Richard Karpinski at the
  43 same address as Thos Sumner.  This includes sample outputs,
  44 bug reports, and criticisms.
  45
  46 You may copy this program freely if you acknowledge its source.
  47 Comments on the Pascal version to NPL, please.
  48
  49 The following is from the introductory commentary from Wichmann's work:
  50
  51 The BASIC program of Kahan is written in Microsoft BASIC using many
  52 facilities which have no exact analogy in Pascal.  The Pascal
  53 version below cannot therefore be exactly the same.  Rather than be
  54 a minimal transcription of the BASIC program, the Pascal coding
  55 follows the conventional style of block-structured languages.  Hence
  56 the Pascal version could be useful in producing versions in other
  57 structured languages.
  58
  59 Rather than use identifiers of minimal length (which therefore have
  60 little mnemonic significance), the Pascal version uses meaningful
  61 identifiers as follows [Note: A few changes have been made for C]:
  62
  63
  64 BASIC   C               BASIC   C               BASIC   C
  65
  66 A                       J                       S    StickyBit
  67 A1   AInverse           J0   NoErrors           T
  68 B    Radix                    [Failure]         T0   Underflow
  69 B1   BInverse           J1   NoErrors           T2   ThirtyTwo
  70 B2   RadixD2                  [SeriousDefect]   T5   OneAndHalf
  71 B9   BMinusU2           J2   NoErrors           T7   TwentySeven
  72 C                             [Defect]          T8   TwoForty
  73 C1   CInverse           J3   NoErrors           U    OneUlp
  74 D                             [Flaw]            U0   UnderflowThreshold
  75 D4   FourD              K    PageNo             U1
  76 E0                      L    Milestone          U2
  77 E1                      M                       V
  78 E2   Exp2               N                       V0
  79 E3                      N1                      V8
  80 E5   MinSqEr            O    Zero               V9
  81 E6   SqEr               O1   One                W
  82 E7   MaxSqEr            O2   Two                X
  83 E8                      O3   Three              X1
  84 E9                      O4   Four               X8
  85 F1   MinusOne           O5   Five               X9   Random1
  86 F2   Half               O8   Eight              Y
  87 F3   Third              O9   Nine               Y1
  88 F6                      P    Precision          Y2
  89 F9                      Q                       Y9   Random2
  90 G1   GMult              Q8                      Z
  91 G2   GDiv               Q9                      Z0   PseudoZero
  92 G3   GAddSub            R                       Z1
  93 H                       R1   RMult              Z2
  94 H1   HInverse           R2   RDiv               Z9
  95 I                       R3   RAddSub
  96 IO   NoTrials           R4   RSqrt
  97 I3   IEEE               R9   Random9
  98
  99 SqRWrng
 100
 101 All the variables in BASIC are true variables and in consequence,
 102 the program is more difficult to follow since the "constants" must
 103 be determined (the glossary is very helpful).  The Pascal version
 104 uses Real constants, but checks are added to ensure that the values
 105 are correctly converted by the compiler.
 106
 107 The major textual change to the Pascal version apart from the
 108 identifiersis that named procedures are used, inserting parameters
 109 wherehelpful.  New procedures are also introduced.  The
 110 correspondence is as follows:
 111
 112
 113 BASIC       Pascal
 114 lines
 115
 116 90- 140   Pause
 117 170- 250   Instructions
 118 380- 460   Heading
 119 480- 670   Characteristics
 120 690- 870   History
 121 2940-2950   Random
 122 3710-3740   NewD
 123 4040-4080   DoesYequalX
 124 4090-4110   PrintIfNPositive
 125 4640-4850   TestPartialUnderflow
 126
 127 */
 128
 129   /* This version of paranoia has been modified to work with GCC's internal
 130      software floating point emulation library, as a sanity check of same.
 131
 132      I'm doing this in C++ so that I can do operator overloading and not
 133      have to modify so damned much of the existing code.  */
 134
 135   extern "C" {
 136 #include <stdio.h>
 137 #include <stddef.h>
 138 #include <limits.h>
 139 #include <string.h>
 140 #include <stdlib.h>
 141 #include <math.h>
 142 #include <unistd.h>
 143 #include <float.h>
 144
 145     /* This part is made all the more awful because many gcc headers are
 146        not prepared at all to be parsed as C++.  The biggest stickler
 147        here is const structure members.  So we include exactly the pieces
 148        that we need.  */
 149
 150 #define GTY(x)
 151
 152 #include "ansidecl.h"
 153 #include "auto-host.h"
 154 #include "hwint.h"
 155
 156 #undef EXTRA_MODES_FILE
 157
 158     struct rtx_def;
 159     typedef struct rtx_def *rtx;
 160     struct rtvec_def;
 161     typedef struct rtvec_def *rtvec;
 162     union tree_node;
 163     typedef union tree_node *tree;
 164
 165 #define DEFTREECODE(SYM, STRING, TYPE, NARGS)   SYM,
 166     enum tree_code {
 167 #include "tree.def"
 168       LAST_AND_UNUSED_TREE_CODE
 169     };
 170 #undef DEFTREECODE
 171
 172 #define ENUM_BITFIELD(X) enum X
 173 #define class klass
 174
 175 #include "real.h"
 176
 177 #undef class
 178   }
 179
 180 /* We never produce signals from the library.  Thus setjmp need do nothing.  */
 181 #undef setjmp
 182 #define setjmp(x)  (0)
 183
 184 static bool verbose = false;
 185 static int verbose_index = 0;
 186
 187 /* ====================================================================== */
 188 /* The implementation of the abstract floating point class based on gcc's
 189    real.c.  I.e. the object of this exercise.  Templated so that we can
 190    all fp sizes.  */
 191
 192 class real_c_float
 193 {
 194  public:
 195   static const enum machine_mode MODE = SFmode;
 196
 197  private:
 198   static const int external_max = 128 / 32;
 199   static const int internal_max
 200     = (sizeof (REAL_VALUE_TYPE) + sizeof (long) + 1) / sizeof (long);
 201   long image[external_max < internal_max ? internal_max : external_max];
 202
 203   void from_long(long);
 204   void from_str(const char *);
 205   void binop(int code, const real_c_float&);
 206   void unop(int code);
 207   bool cmp(int code, const real_c_float&) const;
 208
 209  public:
 210   real_c_float()
 211     { }
 212   real_c_float(long l)
 213     { from_long(l); }
 214   real_c_float(const char *s)
 215     { from_str(s); }
 216   real_c_float(const real_c_float &b)
 217     { memcpy(image, b.image, sizeof(image)); }
 218
 219   const real_c_float& operator= (long l)
 220     { from_long(l); return *this; }
 221   const real_c_float& operator= (const char *s)
 222     { from_str(s); return *this; }
 223   const real_c_float& operator= (const real_c_float &b)
 224     { memcpy(image, b.image, sizeof(image)); return *this; }
 225
 226   const real_c_float& operator+= (const real_c_float &b)
 227     { binop(PLUS_EXPR, b); return *this; }
 228   const real_c_float& operator-= (const real_c_float &b)
 229     { binop(MINUS_EXPR, b); return *this; }
 230   const real_c_float& operator*= (const real_c_float &b)
 231     { binop(MULT_EXPR, b); return *this; }
 232   const real_c_float& operator/= (const real_c_float &b)
 233     { binop(RDIV_EXPR, b); return *this; }
 234
 235   real_c_float operator- () const
 236     { real_c_float r(*this); r.unop(NEGATE_EXPR); return r; }
 237   real_c_float abs () const
 238     { real_c_float r(*this); r.unop(ABS_EXPR); return r; }
 239
 240   bool operator <  (const real_c_float &b) const { return cmp(LT_EXPR, b); }
 241   bool operator <= (const real_c_float &b) const { return cmp(LE_EXPR, b); }
 242   bool operator == (const real_c_float &b) const { return cmp(EQ_EXPR, b); }
 243   bool operator != (const real_c_float &b) const { return cmp(NE_EXPR, b); }
 244   bool operator >= (const real_c_float &b) const { return cmp(GE_EXPR, b); }
 245   bool operator >  (const real_c_float &b) const { return cmp(GT_EXPR, b); }
 246
 247   const char * str () const;
 248   const char * hex () const;
 249   long integer () const;
 250   int exp () const;
 251   void ldexp (int);
 252 };
 253
 254 void
 255 real_c_float::from_long (long l)
 256 {
 257   REAL_VALUE_TYPE f;
 258
 259   real_from_integer (&f, MODE, l, l < 0 ? -1 : 0, 0);
 260   real_to_target (image, &f, MODE);
 261 }
 262
 263 void
 264 real_c_float::from_str (const char *s)
 265 {
 266   REAL_VALUE_TYPE f;
 267   const char *p = s;
 268
 269   if (*p == '-' || *p == '+')
 270     p++;
 271   if (strcasecmp(p, "inf") == 0)
 272     {
 273       real_inf (&f);
 274       if (*s == '-')
 275         real_arithmetic (&f, NEGATE_EXPR, &f, NULL);
 276     }
 277   else if (strcasecmp(p, "nan") == 0)
 278     real_nan (&f, "", 1, MODE);
 279   else
 280     real_from_string (&f, s);
 281
 282   real_to_target (image, &f, MODE);
 283 }
 284
 285 void
 286 real_c_float::binop (int code, const real_c_float &b)
 287 {
 288   REAL_VALUE_TYPE ai, bi, ri;
 289
 290   real_from_target (&ai, image, MODE);
 291   real_from_target (&bi, b.image, MODE);
 292   real_arithmetic (&ri, code, &ai, &bi);
 293   real_to_target (image, &ri, MODE);
 294
 295   if (verbose)
 296     {
 297       char ab[64], bb[64], rb[64];
 298       const real_format *fmt = real_format_for_mode[MODE - QFmode];
 299       const int digits = (fmt->p * fmt->log2_b + 3) / 4;
 300       char symbol_for_code;
 301
 302       real_from_target (&ri, image, MODE);
 303       real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
 304       real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0);
 305       real_to_hexadecimal (rb, &ri, sizeof(rb), digits, 0);
 306
 307       switch (code)
 308         {
 309         case PLUS_EXPR:
 310           symbol_for_code = '+';
 311           break;
 312         case MINUS_EXPR:
 313           symbol_for_code = '-';
 314           break;
 315         case MULT_EXPR:
 316           symbol_for_code = '*';
 317           break;
 318         case RDIV_EXPR:
 319           symbol_for_code = '/';
 320           break;
 321         default:
 322           abort ();
 323         }
 324
 325       fprintf (stderr, "%6d: %s %c %s = %s\n", verbose_index++,
 326                ab, symbol_for_code, bb, rb);
 327     }
 328 }
 329
 330 void
 331 real_c_float::unop (int code)
 332 {
 333   REAL_VALUE_TYPE ai, ri;
 334
 335   real_from_target (&ai, image, MODE);
 336   real_arithmetic (&ri, code, &ai, NULL);
 337   real_to_target (image, &ri, MODE);
 338
 339   if (verbose)
 340     {
 341       char ab[64], rb[64];
 342       const real_format *fmt = real_format_for_mode[MODE - QFmode];
 343       const int digits = (fmt->p * fmt->log2_b + 3) / 4;
 344       const char *symbol_for_code;
 345
 346       real_from_target (&ri, image, MODE);
 347       real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
 348       real_to_hexadecimal (rb, &ri, sizeof(rb), digits, 0);
 349
 350       switch (code)
 351         {
 352         case NEGATE_EXPR:
 353           symbol_for_code = "-";
 354           break;
 355         case ABS_EXPR:
 356           symbol_for_code = "abs ";
 357           break;
 358         default:
 359           abort ();
 360         }
 361
 362       fprintf (stderr, "%6d: %s%s = %s\n", verbose_index++,
 363                symbol_for_code, ab, rb);
 364     }
 365 }
 366
 367 bool
 368 real_c_float::cmp (int code, const real_c_float &b) const
 369 {
 370   REAL_VALUE_TYPE ai, bi;
 371   bool ret;
 372
 373   real_from_target (&ai, image, MODE);
 374   real_from_target (&bi, b.image, MODE);
 375   ret = real_compare (code, &ai, &bi);
 376
 377   if (verbose)
 378     {
 379       char ab[64], bb[64];
 380       const real_format *fmt = real_format_for_mode[MODE - QFmode];
 381       const int digits = (fmt->p * fmt->log2_b + 3) / 4;
 382       const char *symbol_for_code;
 383
 384       real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
 385       real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0);
 386
 387       switch (code)
 388         {
 389         case LT_EXPR:
 390           symbol_for_code = "<";
 391           break;
 392         case LE_EXPR:
 393           symbol_for_code = "<=";
 394           break;
 395         case EQ_EXPR:
 396           symbol_for_code = "==";
 397           break;
 398         case NE_EXPR:
 399           symbol_for_code = "!=";
 400           break;
 401         case GE_EXPR:
 402           symbol_for_code = ">=";
 403           break;
 404         case GT_EXPR:
 405           symbol_for_code = ">";
 406           break;
 407         default:
 408           abort ();
 409         }
 410
 411       fprintf (stderr, "%6d: %s %s %s = %s\n", verbose_index++,
 412                ab, symbol_for_code, bb, (ret ? "true" : "false"));
 413     }
 414
 415   return ret;
 416 }
 417
 418 const char *
 419 real_c_float::str() const
 420 {
 421   REAL_VALUE_TYPE f;
 422   const real_format *fmt = real_format_for_mode[MODE - QFmode];
 423   const int digits = int(fmt->p * fmt->log2_b * .30102999566398119521 + 1);
 424
 425   real_from_target (&f, image, MODE);
 426   char *buf = new char[digits + 10];
 427   real_to_decimal (buf, &f, digits+10, digits, 0);
 428
 429   return buf;
 430 }
 431
 432 const char *
 433 real_c_float::hex() const
 434 {
 435   REAL_VALUE_TYPE f;
 436   const real_format *fmt = real_format_for_mode[MODE - QFmode];
 437   const int digits = (fmt->p * fmt->log2_b + 3) / 4;
 438
 439   real_from_target (&f, image, MODE);
 440   char *buf = new char[digits + 10];
 441   real_to_hexadecimal (buf, &f, digits+10, digits, 0);
 442
 443   return buf;
 444 }
 445
 446 long
 447 real_c_float::integer() const
 448 {
 449   REAL_VALUE_TYPE f;
 450   real_from_target (&f, image, MODE);
 451   return real_to_integer (&f);
 452 }
 453
 454 int
 455 real_c_float::exp() const
 456 {
 457   REAL_VALUE_TYPE f;
 458   real_from_target (&f, image, MODE);
 459   return real_exponent (&f);
 460 }
 461
 462 void
 463 real_c_float::ldexp (int exp)
 464 {
 465   REAL_VALUE_TYPE ai;
 466
 467   real_from_target (&ai, image, MODE);
 468   real_ldexp (&ai, &ai, exp);
 469   real_to_target (image, &ai, MODE);
 470 }
 471
 472 /* ====================================================================== */
 473 /* An implementation of the abstract floating point class that uses native
 474    arithmetic.  Exists for reference and debugging.  */
 475
 476 template<typename T>
 477 class native_float
 478 {
 479  private:
 480   // Force intermediate results back to memory.
 481   volatile T image;
 482
 483   static T from_str (const char *);
 484   static T do_abs (T);
 485   static T verbose_binop (T, char, T, T);
 486   static T verbose_unop (const char *, T, T);
 487   static bool verbose_cmp (T, const char *, T, bool);
 488
 489  public:
 490   native_float()
 491     { }
 492   native_float(long l)
 493     { image = l; }
 494   native_float(const char *s)
 495     { image = from_str(s); }
 496   native_float(const native_float &b)
 497     { image = b.image; }
 498
 499   const native_float& operator= (long l)
 500     { image = l; return *this; }
 501   const native_float& operator= (const char *s)
 502     { image = from_str(s); return *this; }
 503   const native_float& operator= (const native_float &b)
 504     { image = b.image; return *this; }
 505
 506   const native_float& operator+= (const native_float &b)
 507     {
 508       image = verbose_binop(image, '+', b.image, image + b.image);
 509       return *this;
 510     }
 511   const native_float& operator-= (const native_float &b)
 512     {
 513       image = verbose_binop(image, '-', b.image, image - b.image);
 514       return *this;
 515     }
 516   const native_float& operator*= (const native_float &b)
 517     {
 518       image = verbose_binop(image, '*', b.image, image * b.image);
 519       return *this;
 520     }
 521   const native_float& operator/= (const native_float &b)
 522     {
 523       image = verbose_binop(image, '/', b.image, image / b.image);
 524       return *this;
 525     }
 526
 527   native_float operator- () const
 528     {
 529       native_float r;
 530       r.image = verbose_unop("-", image, -image);
 531       return r;
 532     }
 533   native_float abs () const
 534     {
 535       native_float r;
 536       r.image = verbose_unop("abs ", image, do_abs(image));
 537       return r;
 538     }
 539
 540   bool operator <  (const native_float &b) const
 541     { return verbose_cmp(image, "<", b.image, image <  b.image); }
 542   bool operator <= (const native_float &b) const
 543     { return verbose_cmp(image, "<=", b.image, image <= b.image); }
 544   bool operator == (const native_float &b) const
 545     { return verbose_cmp(image, "==", b.image, image == b.image); }
 546   bool operator != (const native_float &b) const
 547     { return verbose_cmp(image, "!=", b.image, image != b.image); }
 548   bool operator >= (const native_float &b) const
 549     { return verbose_cmp(image, ">=", b.image, image >= b.image); }
 550   bool operator >  (const native_float &b) const
 551     { return verbose_cmp(image, ">", b.image, image > b.image); }
 552
 553   const char * str () const;
 554   const char * hex () const;
 555   long integer () const
 556     { return long(image); }
 557   int exp () const;
 558   void ldexp (int);
 559 };
 560
 561 template<typename T>
 562 inline T
 563 native_float<T>::from_str (const char *s)
 564 {
 565   return strtold (s, NULL);
 566 }
 567
 568 template<>
 569 inline float
 570 native_float<float>::from_str (const char *s)
 571 {
 572   return strtof (s, NULL);
 573 }
 574
 575 template<>
 576 inline double
 577 native_float<double>::from_str (const char *s)
 578 {
 579   return strtod (s, NULL);
 580 }
 581
 582 template<typename T>
 583 inline T
 584 native_float<T>::do_abs (T image)
 585 {
 586   return fabsl (image);
 587 }
 588
 589 template<>
 590 inline float
 591 native_float<float>::do_abs (float image)
 592 {
 593   return fabsf (image);
 594 }
 595
 596 template<>
 597 inline double
 598 native_float<double>::do_abs (double image)
 599 {
 600   return fabs (image);
 601 }
 602
 603 template<typename T>
 604 T
 605 native_float<T>::verbose_binop (T a, char symbol, T b, T r)
 606 {
 607   if (verbose)
 608     {
 609       const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1;
 610 #ifdef NO_LONG_DOUBLE
 611       fprintf (stderr, "%6d: %.*a %c %.*a = %.*a\n", verbose_index++,
 612                digits, (double)a, symbol,
 613                digits, (double)b, digits, (double)r);
 614 #else
 615       fprintf (stderr, "%6d: %.*La %c %.*La = %.*La\n", verbose_index++,
 616                digits, (long double)a, symbol,
 617                digits, (long double)b, digits, (long double)r);
 618 #endif
 619     }
 620   return r;
 621 }
 622
 623 template<typename T>
 624 T
 625 native_float<T>::verbose_unop (const char *symbol, T a, T r)
 626 {
 627   if (verbose)
 628     {
 629       const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1;
 630 #ifdef NO_LONG_DOUBLE
 631       fprintf (stderr, "%6d: %s%.*a = %.*a\n", verbose_index++,
 632                symbol, digits, (double)a, digits, (double)r);
 633 #else
 634       fprintf (stderr, "%6d: %s%.*La = %.*La\n", verbose_index++,
 635                symbol, digits, (long double)a, digits, (long double)r);
 636 #endif
 637     }
 638   return r;
 639 }
 640
 641 template<typename T>
 642 bool
 643 native_float<T>::verbose_cmp (T a, const char *symbol, T b, bool r)
 644 {
 645   if (verbose)
 646     {
 647       const int digits = int(sizeof(T) * CHAR_BIT / 4) - 1;
 648 #ifndef NO_LONG_DOUBLE
 649       fprintf (stderr, "%6d: %.*a %s %.*a = %s\n", verbose_index++,
 650                digits, (double)a, symbol,
 651                digits, (double)b, (r ? "true" : "false"));
 652 #else
 653       fprintf (stderr, "%6d: %.*La %s %.*La = %s\n", verbose_index++,
 654                digits, (long double)a, symbol,
 655                digits, (long double)b, (r ? "true" : "false"));
 656 #endif
 657     }
 658   return r;
 659 }
 660
 661 template<typename T>
 662 const char *
 663 native_float<T>::str() const
 664 {
 665   char *buf = new char[50];
 666   const int digits = int(sizeof(T) * CHAR_BIT * .30102999566398119521 + 1);
 667 #ifndef NO_LONG_DOUBLE
 668   sprintf (buf, "%.*e", digits - 1, (double) image);
 669 #else
 670   sprintf (buf, "%.*Le", digits - 1, (long double) image);
 671 #endif
 672   return buf;
 673 }
 674
 675 template<typename T>
 676 const char *
 677 native_float<T>::hex() const
 678 {
 679   char *buf = new char[50];
 680   const int digits = int(sizeof(T) * CHAR_BIT / 4);
 681 #ifndef NO_LONG_DOUBLE
 682   sprintf (buf, "%.*a", digits - 1, (double) image);
 683 #else
 684   sprintf (buf, "%.*La", digits - 1, (long double) image);
 685 #endif
 686   return buf;
 687 }
 688
 689 template<typename T>
 690 int
 691 native_float<T>::exp() const
 692 {
 693   int e;
 694   frexp (image, &e);
 695   return e;
 696 }
 697
 698 template<typename T>
 699 void
 700 native_float<T>::ldexp (int exp)
 701 {
 702   image = ldexpl (image, exp);
 703 }
 704
 705 template<>
 706 void
 707 native_float<float>::ldexp (int exp)
 708 {
 709   image = ldexpf (image, exp);
 710 }
 711
 712 template<>
 713 void
 714 native_float<double>::ldexp (int exp)
 715 {
 716   image = ::ldexp (image, exp);
 717 }
 718
 719 /* ====================================================================== */
 720 /* Some libm routines that Paranoia expects to be available.  */
 721
 722 template<typename FLOAT>
 723 inline FLOAT
 724 FABS (const FLOAT &f)
 725 {
 726   return f.abs();
 727 }
 728
 729 template<typename FLOAT, typename RHS>
 730 inline FLOAT
 731 operator+ (const FLOAT &a, const RHS &b)
 732 {
 733   return FLOAT(a) += FLOAT(b);
 734 }
 735
 736 template<typename FLOAT, typename RHS>
 737 inline FLOAT
 738 operator- (const FLOAT &a, const RHS &b)
 739 {
 740   return FLOAT(a) -= FLOAT(b);
 741 }
 742
 743 template<typename FLOAT, typename RHS>
 744 inline FLOAT
 745 operator* (const FLOAT &a, const RHS &b)
 746 {
 747   return FLOAT(a) *= FLOAT(b);
 748 }
 749
 750 template<typename FLOAT, typename RHS>
 751 inline FLOAT
 752 operator/ (const FLOAT &a, const RHS &b)
 753 {
 754   return FLOAT(a) /= FLOAT(b);
 755 }
 756
 757 template<typename FLOAT>
 758 FLOAT
 759 FLOOR (const FLOAT &f)
 760 {
 761   /* ??? This is only correct when F is representable as an integer.  */
 762   long i = f.integer();
 763   FLOAT r;
 764
 765   r = i;
 766   if (i < 0 && f != r)
 767     r = i - 1;
 768
 769   return r;
 770 }
 771
 772 template<typename FLOAT>
 773 FLOAT
 774 SQRT (const FLOAT &f)
 775 {
 776 #if 0
 777   FLOAT zero = long(0);
 778   FLOAT two = 2;
 779   FLOAT one = 1;
 780   FLOAT diff, diff2;
 781   FLOAT z, t;
 782
 783   if (f == zero)
 784     return zero;
 785   if (f < zero)
 786     return zero / zero;
 787   if (f == one)
 788     return f;
 789
 790   z = f;
 791   z.ldexp (-f.exp() / 2);
 792
 793   diff2 = FABS (z * z - f);
 794   if (diff2 > zero)
 795     while (1)
 796       {
 797         t = (f / (two * z)) + (z / two);
 798         diff = FABS (t * t - f);
 799         if (diff >= diff2)
 800           break;
 801         z = t;
 802         diff2 = diff;
 803       }
 804
 805   return z;
 806 #elif defined(NO_LONG_DOUBLE)
 807   double d;
 808   char buf[64];
 809
 810   d = strtod (f.hex(), NULL);
 811   d = sqrt (d);
 812   sprintf(buf, "%.35a", d);
 813
 814   return FLOAT(buf);
 815 #else
 816   long double ld;
 817   char buf[64];
 818
 819   ld = strtold (f.hex(), NULL);
 820   ld = sqrtl (ld);
 821   sprintf(buf, "%.35La", ld);
 822
 823   return FLOAT(buf);
 824 #endif
 825 }
 826
 827 template<typename FLOAT>
 828 FLOAT
 829 LOG (FLOAT x)
 830 {
 831 #if 0
 832   FLOAT zero = long(0);
 833   FLOAT one = 1;
 834
 835   if (x <= zero)
 836     return zero / zero;
 837   if (x == one)
 838     return zero;
 839
 840   int exp = x.exp() - 1;
 841   x.ldexp(-exp);
 842
 843   FLOAT xm1 = x - one;
 844   FLOAT y = xm1;
 845   long n = 2;
 846
 847   FLOAT sum = xm1;
 848   while (1)
 849     {
 850       y *= xm1;
 851       FLOAT term = y / FLOAT (n);
 852       FLOAT next = sum + term;
 853       if (next == sum)
 854         break;
 855       sum = next;
 856       if (++n == 1000)
 857         break;
 858     }
 859
 860   if (exp)
 861     sum += FLOAT (exp) * FLOAT(".69314718055994530941");
 862
 863   return sum;
 864 #elif defined (NO_LONG_DOUBLE)
 865   double d;
 866   char buf[64];
 867
 868   d = strtod (x.hex(), NULL);
 869   d = log (d);
 870   sprintf(buf, "%.35a", d);
 871
 872   return FLOAT(buf);
 873 #else
 874   long double ld;
 875   char buf[64];
 876
 877   ld = strtold (x.hex(), NULL);
 878   ld = logl (ld);
 879   sprintf(buf, "%.35La", ld);
 880
 881   return FLOAT(buf);
 882 #endif
 883 }
 884
 885 template<typename FLOAT>
 886 FLOAT
 887 EXP (const FLOAT &x)
 888 {
 889   /* Cheat.  */
 890 #ifdef NO_LONG_DOUBLE
 891   double d;
 892   char buf[64];
 893
 894   d = strtod (x.hex(), NULL);
 895   d = exp (d);
 896   sprintf(buf, "%.35a", d);
 897
 898   return FLOAT(buf);
 899 #else
 900   long double ld;
 901   char buf[64];
 902
 903   ld = strtold (x.hex(), NULL);
 904   ld = expl (ld);
 905   sprintf(buf, "%.35La", ld);
 906
 907   return FLOAT(buf);
 908 #endif
 909 }
 910
 911 template<typename FLOAT>
 912 FLOAT
 913 POW (const FLOAT &base, const FLOAT &exp)
 914 {
 915   /* Cheat.  */
 916 #ifdef NO_LONG_DOUBLE
 917   double d1, d2;
 918   char buf[64];
 919
 920   d1 = strtod (base.hex(), NULL);
 921   d2 = strtod (exp.hex(), NULL);
 922   d1 = pow (d1, d2);
 923   sprintf(buf, "%.35a", d1);
 924
 925   return FLOAT(buf);
 926 #else
 927   long double ld1, ld2;
 928   char buf[64];
 929
 930   ld1 = strtold (base.hex(), NULL);
 931   ld2 = strtold (exp.hex(), NULL);
 932   ld1 = powl (ld1, ld2);
 933   sprintf(buf, "%.35La", ld1);
 934
 935   return FLOAT(buf);
 936 #endif
 937 }
 938
 939 /* ====================================================================== */
 940 /* Real Paranoia begins again here.  We wrap the thing in a template so
 941    that we can instantiate it for each floating point type we care for.  */
 942
 943 int NoTrials = 20;              /*Number of tests for commutativity. */
 944 bool do_pause = false;
 945
 946 enum Guard { No, Yes };
 947 enum Rounding { Other, Rounded, Chopped };
 948 enum Class { Failure, Serious, Defect, Flaw };
 949
 950 template<typename FLOAT>
 951 struct Paranoia
 952 {
 953   FLOAT Radix, BInvrse, RadixD2, BMinusU2;
 954
 955   /* Small floating point constants.  */
 956   FLOAT Zero;
 957   FLOAT Half;
 958   FLOAT One;
 959   FLOAT Two;
 960   FLOAT Three;
 961   FLOAT Four;
 962   FLOAT Five;
 963   FLOAT Eight;
 964   FLOAT Nine;
 965   FLOAT TwentySeven;
 966   FLOAT ThirtyTwo;
 967   FLOAT TwoForty;
 968   FLOAT MinusOne;
 969   FLOAT OneAndHalf;
 970
 971   /* Declarations of Variables.  */
 972   int Indx;
 973   char ch[8];
 974   FLOAT AInvrse, A1;
 975   FLOAT C, CInvrse;
 976   FLOAT D, FourD;
 977   FLOAT E0, E1, Exp2, E3, MinSqEr;
 978   FLOAT SqEr, MaxSqEr, E9;
 979   FLOAT Third;
 980   FLOAT F6, F9;
 981   FLOAT H, HInvrse;
 982   int I;
 983   FLOAT StickyBit, J;
 984   FLOAT MyZero;
 985   FLOAT Precision;
 986   FLOAT Q, Q9;
 987   FLOAT R, Random9;
 988   FLOAT T, Underflow, S;
 989   FLOAT OneUlp, UfThold, U1, U2;
 990   FLOAT V, V0, V9;
 991   FLOAT W;
 992   FLOAT X, X1, X2, X8, Random1;
 993   FLOAT Y, Y1, Y2, Random2;
 994   FLOAT Z, PseudoZero, Z1, Z2, Z9;
 995   int ErrCnt[4];
 996   int Milestone;
 997   int PageNo;
 998   int M, N, N1;
 999   Guard GMult, GDiv, GAddSub;
1000   Rounding RMult, RDiv, RAddSub, RSqrt;
1001   int Break, Done, NotMonot, Monot, Anomaly, IEEE, SqRWrng, UfNGrad;
1002
1003   /* Computed constants. */
1004   /*U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
1005   /*U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
1006
1007   int main ();
1008
1009   FLOAT Sign (FLOAT);
1010   FLOAT Random ();
1011   void Pause ();
1012   void BadCond (int, const char *);
1013   void SqXMinX (int);
1014   void TstCond (int, int, const char *);
1015   void notify (const char *);
1016   void IsYeqX ();
1017   void NewD ();
1018   void PrintIfNPositive ();
1019   void SR3750 ();
1020   void TstPtUf ();
1021
1022   // Pretend we're bss.
1023   Paranoia() { memset(this, 0, sizeof (*this)); }
1024 };
1025
1026 template<typename FLOAT>
1027 int
1028 Paranoia<FLOAT>::main()
1029 {
1030   /* First two assignments use integer right-hand sides. */
1031   Zero = long(0);
1032   One = long(1);
1033   Two = long(2);
1034   Three = long(3);
1035   Four = long(4);
1036   Five = long(5);
1037   Eight = long(8);
1038   Nine = long(9);
1039   TwentySeven = long(27);
1040   ThirtyTwo = long(32);
1041   TwoForty = long(240);
1042   MinusOne = long(-1);
1043   Half = "0x1p-1";
1044   OneAndHalf = "0x3p-1";
1045   ErrCnt[Failure] = 0;
1046   ErrCnt[Serious] = 0;
1047   ErrCnt[Defect] = 0;
1048   ErrCnt[Flaw] = 0;
1049   PageNo = 1;
1050         /*=============================================*/
1051   Milestone = 7;
1052         /*=============================================*/
1053   printf ("Program is now RUNNING tests on small integers:\n");
1054
1055   TstCond (Failure, (Zero + Zero == Zero), "0+0 != 0");
1056   TstCond (Failure, (One - One == Zero), "1-1 != 0");
1057   TstCond (Failure, (One > Zero), "1 <= 0");
1058   TstCond (Failure, (One + One == Two), "1+1 != 2");
1059
1060   Z = -Zero;
1061   if (Z != Zero)
1062     {
1063       ErrCnt[Failure] = ErrCnt[Failure] + 1;
1064       printf ("Comparison alleges that -0.0 is Non-zero!\n");
1065       U2 = "0.001";
1066       Radix = 1;
1067       TstPtUf ();
1068     }
1069
1070   TstCond (Failure, (Three == Two + One), "3 != 2+1");
1071   TstCond (Failure, (Four == Three + One), "4 != 3+1");
1072   TstCond (Failure, (Four + Two * (-Two) == Zero), "4 + 2*(-2) != 0");
1073   TstCond (Failure, (Four - Three - One == Zero), "4-3-1 != 0");
1074
1075   TstCond (Failure, (MinusOne == (Zero - One)), "-1 != 0-1");
1076   TstCond (Failure, (MinusOne + One == Zero), "-1+1 != 0");
1077   TstCond (Failure, (One + MinusOne == Zero), "1+(-1) != 0");
1078   TstCond (Failure, (MinusOne + FABS (One) == Zero), "-1+abs(1) != 0");
1079   TstCond (Failure, (MinusOne + MinusOne * MinusOne == Zero),
1080            "-1+(-1)*(-1) != 0");
1081
1082   TstCond (Failure, Half + MinusOne + Half == Zero, "1/2 + (-1) + 1/2 != 0");
1083
1084         /*=============================================*/
1085   Milestone = 10;
1086         /*=============================================*/
1087
1088   TstCond (Failure, (Nine == Three * Three), "9 != 3*3");
1089   TstCond (Failure, (TwentySeven == Nine * Three), "27 != 9*3");
1090   TstCond (Failure, (Eight == Four + Four), "8 != 4+4");
1091   TstCond (Failure, (ThirtyTwo == Eight * Four), "32 != 8*4");
1092   TstCond (Failure, (ThirtyTwo - TwentySeven - Four - One == Zero),
1093            "32-27-4-1 != 0");
1094
1095   TstCond (Failure, Five == Four + One, "5 != 4+1");
1096   TstCond (Failure, TwoForty == Four * Five * Three * Four, "240 != 4*5*3*4");
1097   TstCond (Failure, TwoForty / Three - Four * Four * Five == Zero,
1098            "240/3 - 4*4*5 != 0");
1099   TstCond (Failure, TwoForty / Four - Five * Three * Four == Zero,
1100            "240/4 - 5*3*4 != 0");
1101   TstCond (Failure, TwoForty / Five - Four * Three * Four == Zero,
1102            "240/5 - 4*3*4 != 0");
1103
1104   if (ErrCnt[Failure] == 0)
1105     {
1106       printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
1107       printf ("\n");
1108     }
1109   printf ("Searching for Radix and Precision.\n");
1110   W = One;
1111   do
1112     {
1113       W = W + W;
1114       Y = W + One;
1115       Z = Y - W;
1116       Y = Z - One;
1117     }
1118   while (MinusOne + FABS (Y) < Zero);
1119   /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ... */
1120   Precision = Zero;
1121   Y = One;
1122   do
1123     {
1124       Radix = W + Y;
1125       Y = Y + Y;
1126       Radix = Radix - W;
1127     }
1128   while (Radix == Zero);
1129   if (Radix < Two)
1130     Radix = One;
1131   printf ("Radix = %s .\n", Radix.str());
1132   if (Radix != One)
1133     {
1134       W = One;
1135       do
1136         {
1137           Precision = Precision + One;
1138           W = W * Radix;
1139           Y = W + One;
1140         }
1141       while ((Y - W) == One);
1142     }
1143   /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
1144      ... */
1145   U1 = One / W;
1146   U2 = Radix * U1;
1147   printf ("Closest relative separation found is U1 = %s .\n\n", U1.str());
1148   printf ("Recalculating radix and precision\n ");
1149
1150   /*save old values */
1151   E0 = Radix;
1152   E1 = U1;
1153   E9 = U2;
1154   E3 = Precision;
1155
1156   X = Four / Three;
1157   Third = X - One;
1158   F6 = Half - Third;
1159   X = F6 + F6;
1160   X = FABS (X - Third);
1161   if (X < U2)
1162     X = U2;
1163
1164   /*... now X = (unknown no.) ulps of 1+... */
1165   do
1166     {
1167       U2 = X;
1168       Y = Half * U2 + ThirtyTwo * U2 * U2;
1169       Y = One + Y;
1170       X = Y - One;
1171     }
1172   while (!((U2 <= X) || (X <= Zero)));
1173
1174   /*... now U2 == 1 ulp of 1 + ... */
1175   X = Two / Three;
1176   F6 = X - Half;
1177   Third = F6 + F6;
1178   X = Third - Half;
1179   X = FABS (X + F6);
1180   if (X < U1)
1181     X = U1;
1182
1183   /*... now  X == (unknown no.) ulps of 1 -... */
1184   do
1185     {
1186       U1 = X;
1187       Y = Half * U1 + ThirtyTwo * U1 * U1;
1188       Y = Half - Y;
1189       X = Half + Y;
1190       Y = Half - X;
1191       X = Half + Y;
1192     }
1193   while (!((U1 <= X) || (X <= Zero)));
1194   /*... now U1 == 1 ulp of 1 - ... */
1195   if (U1 == E1)
1196     printf ("confirms closest relative separation U1 .\n");
1197   else
1198     printf ("gets better closest relative separation U1 = %s .\n", U1.str());
1199   W = One / U1;
1200   F9 = (Half - U1) + Half;
1201
1202   Radix = FLOOR (FLOAT ("0.01") + U2 / U1);
1203   if (Radix == E0)
1204     printf ("Radix confirmed.\n");
1205   else
1206     printf ("MYSTERY: recalculated Radix = %s .\n", Radix.str());
1207   TstCond (Defect, Radix <= Eight + Eight,
1208            "Radix is too big: roundoff problems");
1209   TstCond (Flaw, (Radix == Two) || (Radix == 10)
1210            || (Radix == One), "Radix is not as good as 2 or 10");
1211         /*=============================================*/
1212   Milestone = 20;
1213         /*=============================================*/
1214   TstCond (Failure, F9 - Half < Half,
1215            "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
1216   X = F9;
1217   I = 1;
1218   Y = X - Half;
1219   Z = Y - Half;
1220   TstCond (Failure, (X != One)
1221            || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
1222   X = One + U2;
1223   I = 0;
1224         /*=============================================*/
1225   Milestone = 25;
1226         /*=============================================*/
1227   /*... BMinusU2 = nextafter(Radix, 0) */
1228   BMinusU2 = Radix - One;
1229   BMinusU2 = (BMinusU2 - U2) + One;
1230   /* Purify Integers */
1231   if (Radix != One)
1232     {
1233       X = -TwoForty * LOG (U1) / LOG (Radix);
1234       Y = FLOOR (Half + X);
1235       if (FABS (X - Y) * Four < One)
1236         X = Y;
1237       Precision = X / TwoForty;
1238       Y = FLOOR (Half + Precision);
1239       if (FABS (Precision - Y) * TwoForty < Half)
1240         Precision = Y;
1241     }
1242   if ((Precision != FLOOR (Precision)) || (Radix == One))
1243     {
1244       printf ("Precision cannot be characterized by an Integer number\n");
1245       printf
1246         ("of significant digits but, by itself, this is a minor flaw.\n");
1247     }
1248   if (Radix == One)
1249     printf
1250       ("logarithmic encoding has precision characterized solely by U1.\n");
1251   else
1252     printf ("The number of significant digits of the Radix is %s .\n",
1253             Precision.str());
1254   TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
1255            "Precision worse than 5 decimal figures  ");
1256         /*=============================================*/
1257   Milestone = 30;
1258         /*=============================================*/
1259   /* Test for extra-precise subexpressions */
1260   X = FABS (((Four / Three - One) - One / Four) * Three - One / Four);
1261   do
1262     {
1263       Z2 = X;
1264       X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
1265     }
1266   while (!((Z2 <= X) || (X <= Zero)));
1267   X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four);
1268   do
1269     {
1270       Z1 = Z;
1271       Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
1272                         + One / Two)) + One / Two;
1273     }
1274   while (!((Z1 <= Z) || (Z <= Zero)));
1275   do
1276     {
1277       do
1278         {
1279           Y1 = Y;
1280           Y =
1281             (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half)) +
1282             Half;
1283         }
1284       while (!((Y1 <= Y) || (Y <= Zero)));
1285       X1 = X;
1286       X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
1287     }
1288   while (!((X1 <= X) || (X <= Zero)));
1289   if ((X1 != Y1) || (X1 != Z1))
1290     {
1291       BadCond (Serious, "Disagreements among the values X1, Y1, Z1,\n");
1292       printf ("respectively  %s,  %s,  %s,\n", X1.str(), Y1.str(), Z1.str());
1293       printf ("are symptoms of inconsistencies introduced\n");
1294       printf ("by extra-precise evaluation of arithmetic subexpressions.\n");
1295       notify ("Possibly some part of this");
1296       if ((X1 == U1) || (Y1 == U1) || (Z1 == U1))
1297         printf ("That feature is not tested further by this program.\n");
1298     }
1299   else
1300     {
1301       if ((Z1 != U1) || (Z2 != U2))
1302         {
1303           if ((Z1 >= U1) || (Z2 >= U2))
1304             {
1305               BadCond (Failure, "");
1306               notify ("Precision");
1307               printf ("\tU1 = %s, Z1 - U1 = %s\n", U1.str(), (Z1 - U1).str());
1308               printf ("\tU2 = %s, Z2 - U2 = %s\n", U2.str(), (Z2 - U2).str());
1309             }
1310           else
1311             {
1312               if ((Z1 <= Zero) || (Z2 <= Zero))
1313                 {
1314                   printf ("Because of unusual Radix = %s", Radix.str());
1315                   printf (", or exact rational arithmetic a result\n");
1316                   printf ("Z1 = %s, or Z2 = %s ", Z1.str(), Z2.str());
1317                   notify ("of an\nextra-precision");
1318                 }
1319               if (Z1 != Z2 || Z1 > Zero)
1320                 {
1321                   X = Z1 / U1;
1322                   Y = Z2 / U2;
1323                   if (Y > X)
1324                     X = Y;
1325                   Q = -LOG (X);
1326                   printf ("Some subexpressions appear to be calculated "
1327                           "extra precisely\n");
1328                   printf ("with about %s extra B-digits, i.e.\n",
1329                           (Q / LOG (Radix)).str());
1330                   printf ("roughly %s extra significant decimals.\n",
1331                           (Q / LOG (FLOAT (10))).str());
1332                 }
1333               printf
1334                 ("That feature is not tested further by this program.\n");
1335             }
1336         }
1337     }
1338   Pause ();
1339         /*=============================================*/
1340   Milestone = 35;
1341         /*=============================================*/
1342   if (Radix >= Two)
1343     {
1344       X = W / (Radix * Radix);
1345       Y = X + One;
1346       Z = Y - X;
1347       T = Z + U2;
1348       X = T - Z;
1349       TstCond (Failure, X == U2,
1350                "Subtraction is not normalized X=Y,X+Z != Y+Z!");
1351       if (X == U2)
1352         printf ("Subtraction appears to be normalized, as it should be.");
1353     }
1354   printf ("\nChecking for guard digit in *, /, and -.\n");
1355   Y = F9 * One;
1356   Z = One * F9;
1357   X = F9 - Half;
1358   Y = (Y - Half) - X;
1359   Z = (Z - Half) - X;
1360   X = One + U2;
1361   T = X * Radix;
1362   R = Radix * X;
1363   X = T - Radix;
1364   X = X - Radix * U2;
1365   T = R - Radix;
1366   T = T - Radix * U2;
1367   X = X * (Radix - One);
1368   T = T * (Radix - One);
1369   if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero))
1370     GMult = Yes;
1371   else
1372     {
1373       GMult = No;
1374       TstCond (Serious, false, "* lacks a Guard Digit, so 1*X != X");
1375     }
1376   Z = Radix * U2;
1377   X = One + Z;
1378   Y = FABS ((X + Z) - X * X) - U2;
1379   X = One - U2;
1380   Z = FABS ((X - U2) - X * X) - U1;
1381   TstCond (Failure, (Y <= Zero)
1382            && (Z <= Zero), "* gets too many final digits wrong.\n");
1383   Y = One - U2;
1384   X = One + U2;
1385   Z = One / Y;
1386   Y = Z - X;
1387   X = One / Three;
1388   Z = Three / Nine;
1389   X = X - Z;
1390   T = Nine / TwentySeven;
1391   Z = Z - T;
1392   TstCond (Defect, X == Zero && Y == Zero && Z == Zero,
1393            "Division lacks a Guard Digit, so error can exceed 1 ulp\n"
1394            "or  1/3  and  3/9  and  9/27 may disagree");
1395   Y = F9 / One;
1396   X = F9 - Half;
1397   Y = (Y - Half) - X;
1398   X = One + U2;
1399   T = X / One;
1400   X = T - X;
1401   if ((X == Zero) && (Y == Zero) && (Z == Zero))
1402     GDiv = Yes;
1403   else
1404     {
1405       GDiv = No;
1406       TstCond (Serious, false, "Division lacks a Guard Digit, so X/1 != X");
1407     }
1408   X = One / (One + U2);
1409   Y = X - Half - Half;
1410   TstCond (Serious, Y < Zero, "Computed value of 1/1.000..1 >= 1");
1411   X = One - U2;
1412   Y = One + Radix * U2;
1413   Z = X * Radix;
1414   T = Y * Radix;
1415   R = Z / Radix;
1416   StickyBit = T / Radix;
1417   X = R - X;
1418   Y = StickyBit - Y;
1419   TstCond (Failure, X == Zero && Y == Zero,
1420            "* and/or / gets too many last digits wrong");
1421   Y = One - U1;
1422   X = One - F9;
1423   Y = One - Y;
1424   T = Radix - U2;
1425   Z = Radix - BMinusU2;
1426   T = Radix - T;
1427   if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2))
1428     GAddSub = Yes;
1429   else
1430     {
1431       GAddSub = No;
1432       TstCond (Serious, false,
1433                "- lacks Guard Digit, so cancellation is obscured");
1434     }
1435   if (F9 != One && F9 - One >= Zero)
1436     {
1437       BadCond (Serious, "comparison alleges  (1-U1) < 1  although\n");
1438       printf ("  subtraction yields  (1-U1) - 1 = 0 , thereby vitiating\n");
1439       printf ("  such precautions against division by zero as\n");
1440       printf ("  ...  if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
1441     }
1442   if (GMult == Yes && GDiv == Yes && GAddSub == Yes)
1443     printf
1444       ("     *, /, and - appear to have guard digits, as they should.\n");
1445         /*=============================================*/
1446   Milestone = 40;
1447         /*=============================================*/
1448   Pause ();
1449   printf ("Checking rounding on multiply, divide and add/subtract.\n");
1450   RMult = Other;
1451   RDiv = Other;
1452   RAddSub = Other;
1453   RadixD2 = Radix / Two;
1454   A1 = Two;
1455   Done = false;
1456   do
1457     {
1458       AInvrse = Radix;
1459       do
1460         {
1461           X = AInvrse;
1462           AInvrse = AInvrse / A1;
1463         }
1464       while (!(FLOOR (AInvrse) != AInvrse));
1465       Done = (X == One) || (A1 > Three);
1466       if (!Done)
1467         A1 = Nine + One;
1468     }
1469   while (!(Done));
1470   if (X == One)
1471     A1 = Radix;
1472   AInvrse = One / A1;
1473   X = A1;
1474   Y = AInvrse;
1475   Done = false;
1476   do
1477     {
1478       Z = X * Y - Half;
1479       TstCond (Failure, Z == Half, "X * (1/X) differs from 1");
1480       Done = X == Radix;
1481       X = Radix;
1482       Y = One / X;
1483     }
1484   while (!(Done));
1485   Y2 = One + U2;
1486   Y1 = One - U2;
1487   X = OneAndHalf - U2;
1488   Y = OneAndHalf + U2;
1489   Z = (X - U2) * Y2;
1490   T = Y * Y1;
1491   Z = Z - X;
1492   T = T - X;
1493   X = X * Y2;
1494   Y = (Y + U2) * Y1;
1495   X = X - OneAndHalf;
1496   Y = Y - OneAndHalf;
1497   if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero))
1498     {
1499       X = (OneAndHalf + U2) * Y2;
1500       Y = OneAndHalf - U2 - U2;
1501       Z = OneAndHalf + U2 + U2;
1502       T = (OneAndHalf - U2) * Y1;
1503       X = X - (Z + U2);
1504       StickyBit = Y * Y1;
1505       S = Z * Y2;
1506       T = T - Y;
1507       Y = (U2 - Y) + StickyBit;
1508       Z = S - (Z + U2 + U2);
1509       StickyBit = (Y2 + U2) * Y1;
1510       Y1 = Y2 * Y1;
1511       StickyBit = StickyBit - Y2;
1512       Y1 = Y1 - Half;
1513       if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
1514           && (StickyBit == Zero) && (Y1 == Half))
1515         {
1516           RMult = Rounded;
1517           printf ("Multiplication appears to round correctly.\n");
1518         }
1519       else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
1520                && (T < Zero) && (StickyBit + U2 == Zero) && (Y1 < Half))
1521         {
1522           RMult = Chopped;
1523           printf ("Multiplication appears to chop.\n");
1524         }
1525       else
1526         printf ("* is neither chopped nor correctly rounded.\n");
1527       if ((RMult == Rounded) && (GMult == No))
1528         notify ("Multiplication");
1529     }
1530   else
1531     printf ("* is neither chopped nor correctly rounded.\n");
1532         /*=============================================*/
1533   Milestone = 45;
1534         /*=============================================*/
1535   Y2 = One + U2;
1536   Y1 = One - U2;
1537   Z = OneAndHalf + U2 + U2;
1538   X = Z / Y2;
1539   T = OneAndHalf - U2 - U2;
1540   Y = (T - U2) / Y1;
1541   Z = (Z + U2) / Y2;
1542   X = X - OneAndHalf;
1543   Y = Y - T;
1544   T = T / Y1;
1545   Z = Z - (OneAndHalf + U2);
1546   T = (U2 - OneAndHalf) + T;
1547   if (!((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero)))
1548     {
1549       X = OneAndHalf / Y2;
1550       Y = OneAndHalf - U2;
1551       Z = OneAndHalf + U2;
1552       X = X - Y;
1553       T = OneAndHalf / Y1;
1554       Y = Y / Y1;
1555       T = T - (Z + U2);
1556       Y = Y - Z;
1557       Z = Z / Y2;
1558       Y1 = (Y2 + U2) / Y2;
1559       Z = Z - OneAndHalf;
1560       Y2 = Y1 - Y2;
1561       Y1 = (F9 - U1) / F9;
1562       if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
1563           && (Y2 == Zero) && (Y2 == Zero) && (Y1 - Half == F9 - Half))
1564         {
1565           RDiv = Rounded;
1566           printf ("Division appears to round correctly.\n");
1567           if (GDiv == No)
1568             notify ("Division");
1569         }
1570       else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
1571                && (Y2 < Zero) && (Y1 - Half < F9 - Half))
1572         {
1573           RDiv = Chopped;
1574           printf ("Division appears to chop.\n");
1575         }
1576     }
1577   if (RDiv == Other)
1578     printf ("/ is neither chopped nor correctly rounded.\n");
1579   BInvrse = One / Radix;
1580   TstCond (Failure, (BInvrse * Radix - Half == Half),
1581            "Radix * ( 1 / Radix ) differs from 1");
1582         /*=============================================*/
1583   Milestone = 50;
1584         /*=============================================*/
1585   TstCond (Failure, ((F9 + U1) - Half == Half)
1586            && ((BMinusU2 + U2) - One == Radix - One),
1587            "Incomplete carry-propagation in Addition");
1588   X = One - U1 * U1;
1589   Y = One + U2 * (One - U2);
1590   Z = F9 - Half;
1591   X = (X - Half) - Z;
1592   Y = Y - One;
1593   if ((X == Zero) && (Y == Zero))
1594     {
1595       RAddSub = Chopped;
1596       printf ("Add/Subtract appears to be chopped.\n");
1597     }
1598   if (GAddSub == Yes)
1599     {
1600       X = (Half + U2) * U2;
1601       Y = (Half - U2) * U2;
1602       X = One + X;
1603       Y = One + Y;
1604       X = (One + U2) - X;
1605       Y = One - Y;
1606       if ((X == Zero) && (Y == Zero))
1607         {
1608           X = (Half + U2) * U1;
1609           Y = (Half - U2) * U1;
1610           X = One - X;
1611           Y = One - Y;
1612           X = F9 - X;
1613           Y = One - Y;
1614           if ((X == Zero) && (Y == Zero))
1615             {
1616               RAddSub = Rounded;
1617               printf ("Addition/Subtraction appears to round correctly.\n");
1618               if (GAddSub == No)
1619                 notify ("Add/Subtract");
1620             }
1621           else
1622             printf ("Addition/Subtraction neither rounds nor chops.\n");
1623         }
1624       else
1625         printf ("Addition/Subtraction neither rounds nor chops.\n");
1626     }
1627   else
1628     printf ("Addition/Subtraction neither rounds nor chops.\n");
1629   S = One;
1630   X = One + Half * (One + Half);
1631   Y = (One + U2) * Half;
1632   Z = X - Y;
1633   T = Y - X;
1634   StickyBit = Z + T;
1635   if (StickyBit != Zero)
1636     {
1637       S = Zero;
1638       BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n");
1639     }
1640   StickyBit = Zero;
1641   if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
1642       && (RMult == Rounded) && (RDiv == Rounded)
1643       && (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2))
1644     {
1645       printf ("Checking for sticky bit.\n");
1646       X = (Half + U1) * U2;
1647       Y = Half * U2;
1648       Z = One + Y;
1649       T = One + X;
1650       if ((Z - One <= Zero) && (T - One >= U2))
1651         {
1652           Z = T + Y;
1653           Y = Z - X;
1654           if ((Z - T >= U2) && (Y - T == Zero))
1655             {
1656               X = (Half + U1) * U1;
1657               Y = Half * U1;
1658               Z = One - Y;
1659               T = One - X;
1660               if ((Z - One == Zero) && (T - F9 == Zero))
1661                 {
1662                   Z = (Half - U1) * U1;
1663                   T = F9 - Z;
1664                   Q = F9 - Y;
1665                   if ((T - F9 == Zero) && (F9 - U1 - Q == Zero))
1666                     {
1667                       Z = (One + U2) * OneAndHalf;
1668                       T = (OneAndHalf + U2) - Z + U2;
1669                       X = One + Half / Radix;
1670                       Y = One + Radix * U2;
1671                       Z = X * Y;
1672                       if (T == Zero && X + Radix * U2 - Z == Zero)
1673                         {
1674                           if (Radix != Two)
1675                             {
1676                               X = Two + U2;
1677                               Y = X / Two;
1678                               if ((Y - One == Zero))
1679                                 StickyBit = S;
1680                             }
1681                           else
1682                             StickyBit = S;
1683                         }
1684                     }
1685                 }
1686             }
1687         }
1688     }
1689   if (StickyBit == One)
1690     printf ("Sticky bit apparently used correctly.\n");
1691   else
1692     printf ("Sticky bit used incorrectly or not at all.\n");
1693   TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
1694                    RMult == Other || RDiv == Other || RAddSub == Other),
1695            "lack(s) of guard digits or failure(s) to correctly round or chop\n\
1696 (noted above) count as one flaw in the final tally below");
1697         /*=============================================*/
1698   Milestone = 60;
1699         /*=============================================*/
1700   printf ("\n");
1701   printf ("Does Multiplication commute?  ");
1702   printf ("Testing on %d random pairs.\n", NoTrials);
1703   Random9 = SQRT (FLOAT (3));
1704   Random1 = Third;
1705   I = 1;
1706   do
1707     {
1708       X = Random ();
1709       Y = Random ();
1710       Z9 = Y * X;
1711       Z = X * Y;
1712       Z9 = Z - Z9;
1713       I = I + 1;
1714     }
1715   while (!((I > NoTrials) || (Z9 != Zero)));
1716   if (I == NoTrials)
1717     {
1718       Random1 = One + Half / Three;
1719       Random2 = (U2 + U1) + One;
1720       Z = Random1 * Random2;
1721       Y = Random2 * Random1;
1722       Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
1723                                                        Three) * ((U2 + U1) +
1724                                                                  One);
1725     }
1726   if (!((I == NoTrials) || (Z9 == Zero)))
1727     BadCond (Defect, "X * Y == Y * X trial fails.\n");
1728   else
1729     printf ("     No failures found in %d integer pairs.\n", NoTrials);
1730         /*=============================================*/
1731   Milestone = 70;
1732         /*=============================================*/
1733   printf ("\nRunning test of square root(x).\n");
1734   TstCond (Failure, (Zero == SQRT (Zero))
1735            && (-Zero == SQRT (-Zero))
1736            && (One == SQRT (One)), "Square root of 0.0, -0.0 or 1.0 wrong");
1737   MinSqEr = Zero;
1738   MaxSqEr = Zero;
1739   J = Zero;
1740   X = Radix;
1741   OneUlp = U2;
1742   SqXMinX (Serious);
1743   X = BInvrse;
1744   OneUlp = BInvrse * U1;
1745   SqXMinX (Serious);
1746   X = U1;
1747   OneUlp = U1 * U1;
1748   SqXMinX (Serious);
1749   if (J != Zero)
1750     Pause ();
1751   printf ("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
1752   J = Zero;
1753   X = Two;
1754   Y = Radix;
1755   if ((Radix != One))
1756     do
1757       {
1758         X = Y;
1759         Y = Radix * Y;
1760       }
1761     while (!((Y - X >= NoTrials)));
1762   OneUlp = X * U2;
1763   I = 1;
1764   while (I <= NoTrials)
1765     {
1766       X = X + One;
1767       SqXMinX (Defect);
1768       if (J > Zero)
1769         break;
1770       I = I + 1;
1771     }
1772   printf ("Test for sqrt monotonicity.\n");
1773   I = -1;
1774   X = BMinusU2;
1775   Y = Radix;
1776   Z = Radix + Radix * U2;
1777   NotMonot = false;
1778   Monot = false;
1779   while (!(NotMonot || Monot))
1780     {
1781       I = I + 1;
1782       X = SQRT (X);
1783       Q = SQRT (Y);
1784       Z = SQRT (Z);
1785       if ((X > Q) || (Q > Z))
1786         NotMonot = true;
1787       else
1788         {
1789           Q = FLOOR (Q + Half);
1790           if (!(I > 0 || Radix == Q * Q))
1791             Monot = true;
1792           else if (I > 0)
1793             {
1794               if (I > 1)
1795                 Monot = true;
1796               else
1797                 {
1798                   Y = Y * BInvrse;
1799                   X = Y - U1;
1800                   Z = Y + U1;
1801                 }
1802             }
1803           else
1804             {
1805               Y = Q;
1806               X = Y - U2;
1807               Z = Y + U2;
1808             }
1809         }
1810     }
1811   if (Monot)
1812     printf ("sqrt has passed a test for Monotonicity.\n");
1813   else
1814     {
1815       BadCond (Defect, "");
1816       printf ("sqrt(X) is non-monotonic for X near %s .\n", Y.str());
1817     }
1818         /*=============================================*/
1819   Milestone = 110;
1820         /*=============================================*/
1821   printf ("Seeking Underflow thresholds UfThold and E0.\n");
1822   D = U1;
1823   if (Precision != FLOOR (Precision))
1824     {
1825       D = BInvrse;
1826       X = Precision;
1827       do
1828         {
1829           D = D * BInvrse;
1830           X = X - One;
1831         }
1832       while (X > Zero);
1833     }
1834   Y = One;
1835   Z = D;
1836   /* ... D is power of 1/Radix < 1. */
1837   do
1838     {
1839       C = Y;
1840       Y = Z;
1841       Z = Y * Y;
1842     }
1843   while ((Y > Z) && (Z + Z > Z));
1844   Y = C;
1845   Z = Y * D;
1846   do
1847     {
1848       C = Y;
1849       Y = Z;
1850       Z = Y * D;
1851     }
1852   while ((Y > Z) && (Z + Z > Z));
1853   if (Radix < Two)
1854     HInvrse = Two;
1855   else
1856     HInvrse = Radix;
1857   H = One / HInvrse;
1858   /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1859   CInvrse = One / C;
1860   E0 = C;
1861   Z = E0 * H;
1862   /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1863   do
1864     {
1865       Y = E0;
1866       E0 = Z;
1867       Z = E0 * H;
1868     }
1869   while ((E0 > Z) && (Z + Z > Z));
1870   UfThold = E0;
1871   E1 = Zero;
1872   Q = Zero;
1873   E9 = U2;
1874   S = One + E9;
1875   D = C * S;
1876   if (D <= C)
1877     {
1878       E9 = Radix * U2;
1879       S = One + E9;
1880       D = C * S;
1881       if (D <= C)
1882         {
1883           BadCond (Failure,
1884                    "multiplication gets too many last digits wrong.\n");
1885           Underflow = E0;
1886           Y1 = Zero;
1887           PseudoZero = Z;
1888           Pause ();
1889         }
1890     }
1891   else
1892     {
1893       Underflow = D;
1894       PseudoZero = Underflow * H;
1895       UfThold = Zero;
1896       do
1897         {
1898           Y1 = Underflow;
1899           Underflow = PseudoZero;
1900           if (E1 + E1 <= E1)
1901             {
1902               Y2 = Underflow * HInvrse;
1903               E1 = FABS (Y1 - Y2);
1904               Q = Y1;
1905               if ((UfThold == Zero) && (Y1 != Y2))
1906                 UfThold = Y1;
1907             }
1908           PseudoZero = PseudoZero * H;
1909         }
1910       while ((Underflow > PseudoZero)
1911              && (PseudoZero + PseudoZero > PseudoZero));
1912     }
1913   /* Comment line 4530 .. 4560 */
1914   if (PseudoZero != Zero)
1915     {
1916       printf ("\n");
1917       Z = PseudoZero;
1918       /* ... Test PseudoZero for "phoney- zero" violates */
1919       /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1920          ... */
1921       if (PseudoZero <= Zero)
1922         {
1923           BadCond (Failure, "Positive expressions can underflow to an\n");
1924           printf ("allegedly negative value\n");
1925           printf ("PseudoZero that prints out as: %s .\n", PseudoZero.str());
1926           X = -PseudoZero;
1927           if (X <= Zero)
1928             {
1929               printf ("But -PseudoZero, which should be\n");
1930               printf ("positive, isn't; it prints out as  %s .\n", X.str());
1931             }
1932         }
1933       else
1934         {
1935           BadCond (Flaw, "Underflow can stick at an allegedly positive\n");
1936           printf ("value PseudoZero that prints out as %s .\n",
1937                   PseudoZero.str());
1938         }
1939       TstPtUf ();
1940     }
1941         /*=============================================*/
1942   Milestone = 120;
1943         /*=============================================*/
1944   if (CInvrse * Y > CInvrse * Y1)
1945     {
1946       S = H * S;
1947       E0 = Underflow;
1948     }
1949   if (!((E1 == Zero) || (E1 == E0)))
1950     {
1951       BadCond (Defect, "");
1952       if (E1 < E0)
1953         {
1954           printf ("Products underflow at a higher");
1955           printf (" threshold than differences.\n");
1956           if (PseudoZero == Zero)
1957             E0 = E1;
1958         }
1959       else
1960         {
1961           printf ("Difference underflows at a higher");
1962           printf (" threshold than products.\n");
1963         }
1964     }
1965   printf ("Smallest strictly positive number found is E0 = %s .\n", E0.str());
1966   Z = E0;
1967   TstPtUf ();
1968   Underflow = E0;
1969   if (N == 1)
1970     Underflow = Y;
1971   I = 4;
1972   if (E1 == Zero)
1973     I = 3;
1974   if (UfThold == Zero)
1975     I = I - 2;
1976   UfNGrad = true;
1977   switch (I)
1978     {
1979     case 1:
1980       UfThold = Underflow;
1981       if ((CInvrse * Q) != ((CInvrse * Y) * S))
1982         {
1983           UfThold = Y;
1984           BadCond (Failure, "Either accuracy deteriorates as numbers\n");
1985           printf ("approach a threshold = %s\n", UfThold.str());
1986           printf (" coming down from %s\n", C.str());
1987           printf
1988             (" or else multiplication gets too many last digits wrong.\n");
1989         }
1990       Pause ();
1991       break;
1992
1993     case 2:
1994       BadCond (Failure,
1995                "Underflow confuses Comparison, which alleges that\n");
1996       printf ("Q == Y while denying that |Q - Y| == 0; these values\n");
1997       printf ("print out as Q = %s, Y = %s .\n", Q.str(), Y2.str());
1998       printf ("|Q - Y| = %s .\n", FABS (Q - Y2).str());
1999       UfThold = Q;
2000       break;
2001
2002     case 3:
2003       X = X;
2004       break;
2005
2006     case 4:
2007       if ((Q == UfThold) && (E1 == E0) && (FABS (UfThold - E1 / E9) <= E1))
2008         {
2009           UfNGrad = false;
2010           printf ("Underflow is gradual; it incurs Absolute Error =\n");
2011           printf ("(roundoff in UfThold) < E0.\n");
2012           Y = E0 * CInvrse;
2013           Y = Y * (OneAndHalf + U2);
2014           X = CInvrse * (One + U2);
2015           Y = Y / X;
2016           IEEE = (Y == E0);
2017         }
2018     }
2019   if (UfNGrad)
2020     {
2021       printf ("\n");
2022       if (setjmp (ovfl_buf))
2023         {
2024           printf ("Underflow / UfThold failed!\n");
2025           R = H + H;
2026         }
2027       else
2028         R = SQRT (Underflow / UfThold);
2029       if (R <= H)
2030         {
2031           Z = R * UfThold;
2032           X = Z * (One + R * H * (One + H));
2033         }
2034       else
2035         {
2036           Z = UfThold;
2037           X = Z * (One + H * H * (One + H));
2038         }
2039       if (!((X == Z) || (X - Z != Zero)))
2040         {
2041           BadCond (Flaw, "");
2042           printf ("X = %s\n\tis not equal to Z = %s .\n", X.str(), Z.str());
2043           Z9 = X - Z;
2044           printf ("yet X - Z yields %s .\n", Z9.str());
2045           printf ("    Should this NOT signal Underflow, ");
2046           printf ("this is a SERIOUS DEFECT\nthat causes ");
2047           printf ("confusion when innocent statements like\n");;
2048           printf ("    if (X == Z)  ...  else");
2049           printf ("  ... (f(X) - f(Z)) / (X - Z) ...\n");
2050           printf ("encounter Division by Zero although actually\n");
2051           if (setjmp (ovfl_buf))
2052             printf ("X / Z fails!\n");
2053           else
2054             printf ("X / Z = 1 + %s .\n", ((X / Z - Half) - Half).str());
2055         }
2056     }
2057   printf ("The Underflow threshold is %s, below which\n", UfThold.str());
2058   printf ("calculation may suffer larger Relative error than ");
2059   printf ("merely roundoff.\n");
2060   Y2 = U1 * U1;
2061   Y = Y2 * Y2;
2062   Y2 = Y * U1;
2063   if (Y2 <= UfThold)
2064     {
2065       if (Y > E0)
2066         {
2067           BadCond (Defect, "");
2068           I = 5;
2069         }
2070       else
2071         {
2072           BadCond (Serious, "");
2073           I = 4;
2074         }
2075       printf ("Range is too narrow; U1^%d Underflows.\n", I);
2076     }
2077         /*=============================================*/
2078   Milestone = 130;
2079         /*=============================================*/
2080   Y = -FLOOR (Half - TwoForty * LOG (UfThold) / LOG (HInvrse)) / TwoForty;
2081   Y2 = Y + Y;
2082   printf ("Since underflow occurs below the threshold\n");
2083   printf ("UfThold = (%s) ^ (%s)\nonly underflow ", HInvrse.str(), Y.str());
2084   printf ("should afflict the expression\n\t(%s) ^ (%s);\n",
2085           HInvrse.str(), Y2.str());
2086   printf ("actually calculating yields:");
2087   if (setjmp (ovfl_buf))
2088     {
2089       BadCond (Serious, "trap on underflow.\n");
2090     }
2091   else
2092     {
2093       V9 = POW (HInvrse, Y2);
2094       printf (" %s .\n", V9.str());
2095       if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold)))
2096         {
2097           BadCond (Serious, "this is not between 0 and underflow\n");
2098           printf ("   threshold = %s .\n", UfThold.str());
2099         }
2100       else if (!(V9 > UfThold * (One + E9)))
2101         printf ("This computed value is O.K.\n");
2102       else
2103         {
2104           BadCond (Defect, "this is not between 0 and underflow\n");
2105           printf ("   threshold = %s .\n", UfThold.str());
2106         }
2107     }
2108         /*=============================================*/
2109   Milestone = 160;
2110         /*=============================================*/
2111   Pause ();
2112   printf ("Searching for Overflow threshold:\n");
2113   printf ("This may generate an error.\n");
2114   Y = -CInvrse;
2115   V9 = HInvrse * Y;
2116   if (setjmp (ovfl_buf))
2117     {
2118       I = 0;
2119       V9 = Y;
2120       goto overflow;
2121     }
2122   do
2123     {
2124       V = Y;
2125       Y = V9;
2126       V9 = HInvrse * Y;
2127     }
2128   while (V9 < Y);
2129   I = 1;
2130 overflow:
2131   Z = V9;
2132   printf ("Can `Z = -Y' overflow?\n");
2133   printf ("Trying it on Y = %s .\n", Y.str());
2134   V9 = -Y;
2135   V0 = V9;
2136   if (V - Y == V + V0)
2137     printf ("Seems O.K.\n");
2138   else
2139     {
2140       printf ("finds a ");
2141       BadCond (Flaw, "-(-Y) differs from Y.\n");
2142     }
2143   if (Z != Y)
2144     {
2145       BadCond (Serious, "");
2146       printf ("overflow past %s\n\tshrinks to %s .\n", Y.str(), Z.str());
2147     }
2148   if (I)
2149     {
2150       Y = V * (HInvrse * U2 - HInvrse);
2151       Z = Y + ((One - HInvrse) * U2) * V;
2152       if (Z < V0)
2153         Y = Z;
2154       if (Y < V0)
2155         V = Y;
2156       if (V0 - V < V0)
2157         V = V0;
2158     }
2159   else
2160     {
2161       V = Y * (HInvrse * U2 - HInvrse);
2162       V = V + ((One - HInvrse) * U2) * Y;
2163     }
2164   printf ("Overflow threshold is V  = %s .\n", V.str());
2165   if (I)
2166     printf ("Overflow saturates at V0 = %s .\n", V0.str());
2167   else
2168     printf ("There is no saturation value because "
2169             "the system traps on overflow.\n");
2170   V9 = V * One;
2171   printf ("No Overflow should be signaled for V * 1 = %s\n", V9.str());
2172   V9 = V / One;
2173   printf ("                           nor for V / 1 = %s.\n", V9.str());
2174   printf ("Any overflow signal separating this * from the one\n");
2175   printf ("above is a DEFECT.\n");
2176         /*=============================================*/
2177   Milestone = 170;
2178         /*=============================================*/
2179   if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V))
2180     {
2181       BadCond (Failure, "Comparisons involving ");
2182       printf ("+-%s, +-%s\nand +-%s are confused by Overflow.",
2183               V.str(), V0.str(), UfThold.str());
2184     }
2185         /*=============================================*/
2186   Milestone = 175;
2187         /*=============================================*/
2188   printf ("\n");
2189   for (Indx = 1; Indx <= 3; ++Indx)
2190     {
2191       switch (Indx)
2192         {
2193         case 1:
2194           Z = UfThold;
2195           break;
2196         case 2:
2197           Z = E0;
2198           break;
2199         case 3:
2200           Z = PseudoZero;
2201           break;
2202         }
2203       if (Z != Zero)
2204         {
2205           V9 = SQRT (Z);
2206           Y = V9 * V9;
2207           if (Y / (One - Radix * E9) < Z || Y > (One + Radix * E9) * Z)
2208             {                   /* dgh: + E9 --> * E9 */
2209               if (V9 > U1)
2210                 BadCond (Serious, "");
2211               else
2212                 BadCond (Defect, "");
2213               printf ("Comparison alleges that what prints as Z = %s\n",
2214                       Z.str());
2215               printf (" is too far from sqrt(Z) ^ 2 = %s .\n", Y.str());
2216             }
2217         }
2218     }
2219         /*=============================================*/
2220   Milestone = 180;
2221         /*=============================================*/
2222   for (Indx = 1; Indx <= 2; ++Indx)
2223     {
2224       if (Indx == 1)
2225         Z = V;
2226       else
2227         Z = V0;
2228       V9 = SQRT (Z);
2229       X = (One - Radix * E9) * V9;
2230       V9 = V9 * X;
2231       if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z)))
2232         {
2233           Y = V9;
2234           if (X < W)
2235             BadCond (Serious, "");
2236           else
2237             BadCond (Defect, "");
2238           printf ("Comparison alleges that Z = %s\n", Z.str());
2239           printf (" is too far from sqrt(Z) ^ 2 (%s) .\n", Y.str());
2240         }
2241     }
2242         /*=============================================*/
2243   Milestone = 190;
2244         /*=============================================*/
2245   Pause ();
2246   X = UfThold * V;
2247   Y = Radix * Radix;
2248   if (X * Y < One || X > Y)
2249     {
2250       if (X * Y < U1 || X > Y / U1)
2251         BadCond (Defect, "Badly");
2252       else
2253         BadCond (Flaw, "");
2254
2255       printf (" unbalanced range; UfThold * V = %s\n\t%s\n",
2256               X.str(), "is too far from 1.\n");
2257     }
2258         /*=============================================*/
2259   Milestone = 200;
2260         /*=============================================*/
2261   for (Indx = 1; Indx <= 5; ++Indx)
2262     {
2263       X = F9;
2264       switch (Indx)
2265         {
2266         case 2:
2267           X = One + U2;
2268           break;
2269         case 3:
2270           X = V;
2271           break;
2272         case 4:
2273           X = UfThold;
2274           break;
2275         case 5:
2276           X = Radix;
2277         }
2278       Y = X;
2279       if (setjmp (ovfl_buf))
2280         printf ("  X / X  traps when X = %s\n", X.str());
2281       else
2282         {
2283           V9 = (Y / X - Half) - Half;
2284           if (V9 == Zero)
2285             continue;
2286           if (V9 == -U1 && Indx < 5)
2287             BadCond (Flaw, "");
2288           else
2289             BadCond (Serious, "");
2290           printf ("  X / X differs from 1 when X = %s\n", X.str());
2291           printf ("  instead, X / X - 1/2 - 1/2 = %s .\n", V9.str());
2292         }
2293     }
2294         /*=============================================*/
2295   Milestone = 210;
2296         /*=============================================*/
2297   MyZero = Zero;
2298   printf ("\n");
2299   printf ("What message and/or values does Division by Zero produce?\n");
2300   printf ("    Trying to compute 1 / 0 produces ...");
2301   if (!setjmp (ovfl_buf))
2302     printf ("  %s .\n", (One / MyZero).str());
2303   printf ("\n    Trying to compute 0 / 0 produces ...");
2304   if (!setjmp (ovfl_buf))
2305     printf ("  %s .\n", (Zero / MyZero).str());
2306         /*=============================================*/
2307   Milestone = 220;
2308         /*=============================================*/
2309   Pause ();
2310   printf ("\n");
2311   {
2312     static const char *msg[] = {
2313       "FAILUREs  encountered =",
2314       "SERIOUS DEFECTs  discovered =",
2315       "DEFECTs  discovered =",
2316       "FLAWs  discovered ="
2317     };
2318     int i;
2319     for (i = 0; i < 4; i++)
2320       if (ErrCnt[i])
2321         printf ("The number of  %-29s %d.\n", msg[i], ErrCnt[i]);
2322   }
2323   printf ("\n");
2324   if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] + ErrCnt[Flaw]) > 0)
2325     {
2326       if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] == 0)
2327           && (ErrCnt[Flaw] > 0))
2328         {
2329           printf ("The arithmetic diagnosed seems ");
2330           printf ("Satisfactory though flawed.\n");
2331         }
2332       if ((ErrCnt[Failure] + ErrCnt[Serious] == 0) && (ErrCnt[Defect] > 0))
2333         {
2334           printf ("The arithmetic diagnosed may be Acceptable\n");
2335           printf ("despite inconvenient Defects.\n");
2336         }
2337       if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0)
2338         {
2339           printf ("The arithmetic diagnosed has ");
2340           printf ("unacceptable Serious Defects.\n");
2341         }
2342       if (ErrCnt[Failure] > 0)
2343         {
2344           printf ("Potentially fatal FAILURE may have spoiled this");
2345           printf (" program's subsequent diagnoses.\n");
2346         }
2347     }
2348   else
2349     {
2350       printf ("No failures, defects nor flaws have been discovered.\n");
2351       if (!((RMult == Rounded) && (RDiv == Rounded)
2352             && (RAddSub == Rounded) && (RSqrt == Rounded)))
2353         printf ("The arithmetic diagnosed seems Satisfactory.\n");
2354       else
2355         {
2356           if (StickyBit >= One &&
2357               (Radix - Two) * (Radix - Nine - One) == Zero)
2358             {
2359               printf ("Rounding appears to conform to ");
2360               printf ("the proposed IEEE standard P");
2361               if ((Radix == Two) &&
2362                   ((Precision - Four * Three * Two) *
2363                    (Precision - TwentySeven - TwentySeven + One) == Zero))
2364                 printf ("754");
2365               else
2366                 printf ("854");
2367               if (IEEE)
2368                 printf (".\n");
2369               else
2370                 {
2371                   printf (",\nexcept for possibly Double Rounding");
2372                   printf (" during Gradual Underflow.\n");
2373                 }
2374             }
2375           printf ("The arithmetic diagnosed appears to be Excellent!\n");
2376         }
2377     }
2378   printf ("END OF TEST.\n");
2379   return 0;
2380 }
2381
2382 template<typename FLOAT>
2383 FLOAT
2384 Paranoia<FLOAT>::Sign (FLOAT X)
2385 {
2386   return X >= FLOAT (long (0)) ? 1 : -1;
2387 }
2388
2389 template<typename FLOAT>
2390 void
2391 Paranoia<FLOAT>::Pause ()
2392 {
2393   if (do_pause)
2394     {
2395       fputs ("Press return...", stdout);
2396       fflush (stdout);
2397       getchar();
2398     }
2399   printf ("\nDiagnosis resumes after milestone Number %d", Milestone);
2400   printf ("          Page: %d\n\n", PageNo);
2401   ++Milestone;
2402   ++PageNo;
2403 }
2404
2405 template<typename FLOAT>
2406 void
2407 Paranoia<FLOAT>::TstCond (int K, int Valid, const char *T)
2408 {
2409   if (!Valid)
2410     {
2411       BadCond (K, T);
2412       printf (".\n");
2413     }
2414 }
2415
2416 template<typename FLOAT>
2417 void
2418 Paranoia<FLOAT>::BadCond (int K, const char *T)
2419 {
2420   static const char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
2421
2422   ErrCnt[K] = ErrCnt[K] + 1;
2423   printf ("%s:  %s", msg[K], T);
2424 }
2425
2426 /* Random computes
2427      X = (Random1 + Random9)^5
2428      Random1 = X - FLOOR(X) + 0.000005 * X;
2429    and returns the new value of Random1.  */
2430
2431 template<typename FLOAT>
2432 FLOAT
2433 Paranoia<FLOAT>::Random ()
2434 {
2435   FLOAT X, Y;
2436
2437   X = Random1 + Random9;
2438   Y = X * X;
2439   Y = Y * Y;
2440   X = X * Y;
2441   Y = X - FLOOR (X);
2442   Random1 = Y + X * FLOAT ("0.000005");
2443   return (Random1);
2444 }
2445
2446 template<typename FLOAT>
2447 void
2448 Paranoia<FLOAT>::SqXMinX (int ErrKind)
2449 {
2450   FLOAT XA, XB;
2451
2452   XB = X * BInvrse;
2453   XA = X - XB;
2454   SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp;
2455   if (SqEr != Zero)
2456     {
2457       if (SqEr < MinSqEr)
2458         MinSqEr = SqEr;
2459       if (SqEr > MaxSqEr)
2460         MaxSqEr = SqEr;
2461       J = J + 1;
2462       BadCond (ErrKind, "\n");
2463       printf ("sqrt(%s) - %s  = %s\n", (X * X).str(), X.str(),
2464               (OneUlp * SqEr).str());
2465       printf ("\tinstead of correct value 0 .\n");
2466     }
2467 }
2468
2469 template<typename FLOAT>
2470 void
2471 Paranoia<FLOAT>::NewD ()
2472 {
2473   X = Z1 * Q;
2474   X = FLOOR (Half - X / Radix) * Radix + X;
2475   Q = (Q - X * Z) / Radix + X * X * (D / Radix);
2476   Z = Z - Two * X * D;
2477   if (Z <= Zero)
2478     {
2479       Z = -Z;
2480       Z1 = -Z1;
2481     }
2482   D = Radix * D;
2483 }
2484
2485 template<typename FLOAT>
2486 void
2487 Paranoia<FLOAT>::SR3750 ()
2488 {
2489   if (!((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2)))
2490     {
2491       I = I + 1;
2492       X2 = SQRT (X * D);
2493       Y2 = (X2 - Z2) - (Y - Z2);
2494       X2 = X8 / (Y - Half);
2495       X2 = X2 - Half * X2 * X2;
2496       SqEr = (Y2 + Half) + (Half - X2);
2497       if (SqEr < MinSqEr)
2498         MinSqEr = SqEr;
2499       SqEr = Y2 - X2;
2500       if (SqEr > MaxSqEr)
2501         MaxSqEr = SqEr;
2502     }
2503 }
2504
2505 template<typename FLOAT>
2506 void
2507 Paranoia<FLOAT>::IsYeqX ()
2508 {
2509   if (Y != X)
2510     {
2511       if (N <= 0)
2512         {
2513           if (Z == Zero && Q <= Zero)
2514             printf ("WARNING:  computing\n");
2515           else
2516             BadCond (Defect, "computing\n");
2517           printf ("\t(%s) ^ (%s)\n", Z.str(), Q.str());
2518           printf ("\tyielded %s;\n", Y.str());
2519           printf ("\twhich compared unequal to correct %s ;\n", X.str());
2520           printf ("\t\tthey differ by %s .\n", (Y - X).str());
2521         }
2522       N = N + 1;                /* ... count discrepancies. */
2523     }
2524 }
2525
2526 template<typename FLOAT>
2527 void
2528 Paranoia<FLOAT>::PrintIfNPositive ()
2529 {
2530   if (N > 0)
2531     printf ("Similar discrepancies have occurred %d times.\n", N);
2532 }
2533
2534 template<typename FLOAT>
2535 void
2536 Paranoia<FLOAT>::TstPtUf ()
2537 {
2538   N = 0;
2539   if (Z != Zero)
2540     {
2541       printf ("Since comparison denies Z = 0, evaluating ");
2542       printf ("(Z + Z) / Z should be safe.\n");
2543       if (setjmp (ovfl_buf))
2544         goto very_serious;
2545       Q9 = (Z + Z) / Z;
2546       printf ("What the machine gets for (Z + Z) / Z is %s .\n", Q9.str());
2547       if (FABS (Q9 - Two) < Radix * U2)
2548         {
2549           printf ("This is O.K., provided Over/Underflow");
2550           printf (" has NOT just been signaled.\n");
2551         }
2552       else
2553         {
2554           if ((Q9 < One) || (Q9 > Two))
2555             {
2556             very_serious:
2557               N = 1;
2558               ErrCnt[Serious] = ErrCnt[Serious] + 1;
2559               printf ("This is a VERY SERIOUS DEFECT!\n");
2560             }
2561           else
2562             {
2563               N = 1;
2564               ErrCnt[Defect] = ErrCnt[Defect] + 1;
2565               printf ("This is a DEFECT!\n");
2566             }
2567         }
2568       V9 = Z * One;
2569       Random1 = V9;
2570       V9 = One * Z;
2571       Random2 = V9;
2572       V9 = Z / One;
2573       if ((Z == Random1) && (Z == Random2) && (Z == V9))
2574         {
2575           if (N > 0)
2576             Pause ();
2577         }
2578       else
2579         {
2580           N = 1;
2581           BadCond (Defect, "What prints as Z = ");
2582           printf ("%s\n\tcompares different from  ", Z.str());
2583           if (Z != Random1)
2584             printf ("Z * 1 = %s ", Random1.str());
2585           if (!((Z == Random2) || (Random2 == Random1)))
2586             printf ("1 * Z == %s\n", Random2.str());
2587           if (!(Z == V9))
2588             printf ("Z / 1 = %s\n", V9.str());
2589           if (Random2 != Random1)
2590             {
2591               ErrCnt[Defect] = ErrCnt[Defect] + 1;
2592               BadCond (Defect, "Multiplication does not commute!\n");
2593               printf ("\tComparison alleges that 1 * Z = %s\n", Random2.str());
2594               printf ("\tdiffers from Z * 1 = %s\n", Random1.str());
2595             }
2596           Pause ();
2597         }
2598     }
2599 }
2600
2601 template<typename FLOAT>
2602 void
2603 Paranoia<FLOAT>::notify (const char *s)
2604 {
2605   printf ("%s test appears to be inconsistent...\n", s);
2606   printf ("   PLEASE NOTIFY KARPINKSI!\n");
2607 }
2608
2609 /* ====================================================================== */
2610
2611 int main(int ac, char **av)
2612 {
2613   setbuf(stdout, NULL);
2614   setbuf(stderr, NULL);
2615
2616   while (1)
2617     switch (getopt (ac, av, "pvg:fdl"))
2618       {
2619       case -1:
2620         return 0;
2621       case 'p':
2622         do_pause = true;
2623         break;
2624       case 'v':
2625         verbose = true;
2626         break;
2627       case 'g':
2628         {
2629           static const struct {
2630             const char *name;
2631             const struct real_format *fmt;
2632           } fmts[] = {
2633 #define F(x) { #x, &x##_format }
2634             F(ieee_single),
2635             F(ieee_double),
2636             F(ieee_extended_motorola),
2637             F(ieee_extended_intel_96),
2638             F(ieee_extended_intel_128),
2639             F(ibm_extended),
2640             F(ieee_quad),
2641             F(vax_f),
2642             F(vax_d),
2643             F(vax_g),
2644             F(i370_single),
2645             F(i370_double),
2646             F(c4x_single),
2647             F(c4x_extended),
2648             F(real_internal),
2649 #undef F
2650           };
2651
2652           int i, n = sizeof (fmts)/sizeof(*fmts);
2653
2654           for (i = 0; i < n; ++i)
2655             if (strcmp (fmts[i].name, optarg) == 0)
2656               break;
2657
2658           if (i == n)
2659             {
2660               printf ("Unknown implementation \"%s\"; "
2661                       "available implementations:\n", optarg);
2662               for (i = 0; i < n; ++i)
2663                 printf ("\t%s\n", fmts[i].name);
2664               return 1;
2665             }
2666
2667           // We cheat and use the same mode all the time, but vary
2668           // the format used for that mode.
2669           real_format_for_mode[int(real_c_float::MODE) - int(QFmode)]
2670             = fmts[i].fmt;
2671
2672           Paranoia<real_c_float>().main();
2673           break;
2674         }
2675
2676       case 'f':
2677         Paranoia < native_float<float> >().main();
2678         break;
2679       case 'd':
2680         Paranoia < native_float<double> >().main();
2681         break;
2682       case 'l':
2683 #ifndef NO_LONG_DOUBLE
2684         Paranoia < native_float<long double> >().main();
2685 #endif
2686         break;
2687
2688       case '?':
2689         puts ("-p\tpause between pages");
2690         puts ("-g<FMT>\treal.c implementation FMT");
2691         puts ("-f\tnative float");
2692         puts ("-d\tnative double");
2693         puts ("-l\tnative long double");
2694         return 0;
2695       }
2696 }
2697
2698 /* GCC stuff referenced by real.o.  */
2699
2700 extern "C" void
2701 fancy_abort ()
2702 {
2703   abort ();
2704 }
2705
2706 int target_flags = 0;
2707
2708 extern "C" int
2709 floor_log2_wide (unsigned HOST_WIDE_INT x)
2710 {
2711   int log = -1;
2712   while (x != 0)
2713     log++,
2714     x >>= 1;
2715   return log;
2716 }