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