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