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 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 183static bool verbose = false; 184static 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 191class 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 253void 254real_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 262void 263real_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 284void 285real_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 329void 330real_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 366bool 367real_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 417const char * 418real_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 431const char * 432real_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 445long 446real_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 453int 454real_c_float::exp() const 455{ 456 REAL_VALUE_TYPE f; 457 real_from_target (&f, image, MODE); 458 return real_exponent (&f); 459} 460 461void 462real_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 475template<typename T> 476class 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 560template<typename T> 561inline T 562native_float<T>::from_str (const char *s) 563{ 564 return strtold (s, NULL); 565} 566 567template<> 568inline float 569native_float<float>::from_str (const char *s) 570{ 571 return strtof (s, NULL); 572} 573 574template<> 575inline double 576native_float<double>::from_str (const char *s) 577{ 578 return strtod (s, NULL); 579} 580 581template<typename T> 582inline T 583native_float<T>::do_abs (T image) 584{ 585 return fabsl (image); 586} 587 588template<> 589inline float 590native_float<float>::do_abs (float image) 591{ 592 return fabsf (image); 593} 594 595template<> 596inline double 597native_float<double>::do_abs (double image) 598{ 599 return fabs (image); 600} 601 602template<typename T> 603T 604native_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 622template<typename T> 623T 624native_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 640template<typename T> 641bool 642native_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 660template<typename T> 661const char * 662native_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 674template<typename T> 675const char * 676native_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 688template<typename T> 689int 690native_float<T>::exp() const 691{ 692 int e; 693 frexp (image, &e); 694 return e; 695} 696 697template<typename T> 698void 699native_float<T>::ldexp (int exp) 700{ 701 image = ldexpl (image, exp); 702} 703 704template<> 705void 706native_float<float>::ldexp (int exp) 707{ 708 image = ldexpf (image, exp); 709} 710 711template<> 712void 713native_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 721template<typename FLOAT> 722inline FLOAT 723FABS (const FLOAT &f) 724{ 725 return f.abs(); 726} 727 728template<typename FLOAT, typename RHS> 729inline FLOAT 730operator+ (const FLOAT &a, const RHS &b) 731{ 732 return FLOAT(a) += FLOAT(b); 733} 734 735template<typename FLOAT, typename RHS> 736inline FLOAT 737operator- (const FLOAT &a, const RHS &b) 738{ 739 return FLOAT(a) -= FLOAT(b); 740} 741 742template<typename FLOAT, typename RHS> 743inline FLOAT 744operator* (const FLOAT &a, const RHS &b) 745{ 746 return FLOAT(a) *= FLOAT(b); 747} 748 749template<typename FLOAT, typename RHS> 750inline FLOAT 751operator/ (const FLOAT &a, const RHS &b) 752{ 753 return FLOAT(a) /= FLOAT(b); 754} 755 756template<typename FLOAT> 757FLOAT 758FLOOR (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 771template<typename FLOAT> 772FLOAT 773SQRT (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 826template<typename FLOAT> 827FLOAT 828LOG (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 884template<typename FLOAT> 885FLOAT 886EXP (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 910template<typename FLOAT> 911FLOAT 912POW (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 942int NoTrials = 20; /*Number of tests for commutativity. */ 943bool do_pause = false; 944 945enum Guard { No, Yes }; 946enum Rounding { Other, Rounded, Chopped }; 947enum Class { Failure, Serious, Defect, Flaw }; 948 949template<typename FLOAT> 950struct 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 1025template<typename FLOAT> 1026int 1027Paranoia<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; 2129overflow: 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 2381template<typename FLOAT> 2382FLOAT 2383Paranoia<FLOAT>::Sign (FLOAT X) 2384{ 2385 return X >= FLOAT (long (0)) ? 1 : -1; 2386} 2387 2388template<typename FLOAT> 2389void 2390Paranoia<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 2404template<typename FLOAT> 2405void 2406Paranoia<FLOAT>::TstCond (int K, int Valid, const char *T) 2407{ 2408 if (!Valid) 2409 { 2410 BadCond (K, T); 2411 printf (".\n"); 2412 } 2413} 2414 2415template<typename FLOAT> 2416void 2417Paranoia<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 2430template<typename FLOAT> 2431FLOAT 2432Paranoia<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 2445template<typename FLOAT> 2446void 2447Paranoia<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 2468template<typename FLOAT> 2469void 2470Paranoia<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 2484template<typename FLOAT> 2485void 2486Paranoia<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 2504template<typename FLOAT> 2505void 2506Paranoia<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 2525template<typename FLOAT> 2526void 2527Paranoia<FLOAT>::PrintIfNPositive () 2528{ 2529 if (N > 0) 2530 printf ("Similar discrepancies have occurred %d times.\n", N); 2531} 2532 2533template<typename FLOAT> 2534void 2535Paranoia<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 2600template<typename FLOAT> 2601void 2602Paranoia<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 2610int 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 2697extern "C" void 2698fancy_abort () 2699{ 2700 abort (); 2701} 2702 2703int target_flags = 0; 2704 2705extern "C" int 2706floor_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} 2714