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