Mark ChangeLog
[official-gcc.git] / contrib / paranoia.cc
blobad9f8c8fda7864262b6be3481e6c5518fda959d2
1 /* A C version of Kahan's Floating Point Test "Paranoia"
3 Thos Sumner, UCSF, Feb. 1985
4 David Gay, BTL, Jan. 1986
6 This is a rewrite from the Pascal version by
8 B. A. Wichmann, 18 Jan. 1985
10 (and does NOT exhibit good C programming style).
12 Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
14 (C) Apr 19 1983 in BASIC version by:
15 Professor W. M. Kahan,
16 567 Evans Hall
17 Electrical Engineering & Computer Science Dept.
18 University of California
19 Berkeley, California 94720
20 USA
22 converted to Pascal by:
23 B. A. Wichmann
24 National Physical Laboratory
25 Teddington Middx
26 TW11 OLW
29 converted to C by:
31 David M. Gay and Thos Sumner
32 AT&T Bell Labs Computer Center, Rm. U-76
33 600 Mountain Avenue University of California
34 Murray Hill, NJ 07974 San Francisco, CA 94143
35 USA USA
37 with simultaneous corrections to the Pascal source (reflected
38 in the Pascal source available over netlib).
39 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
41 Reports of results on various systems from all the versions
42 of Paranoia are being collected by Richard Karpinski at the
43 same address as Thos Sumner. This includes sample outputs,
44 bug reports, and criticisms.
46 You may copy this program freely if you acknowledge its source.
47 Comments on the Pascal version to NPL, please.
49 The following is from the introductory commentary from Wichmann's work:
51 The BASIC program of Kahan is written in Microsoft BASIC using many
52 facilities which have no exact analogy in Pascal. The Pascal
53 version below cannot therefore be exactly the same. Rather than be
54 a minimal transcription of the BASIC program, the Pascal coding
55 follows the conventional style of block-structured languages. Hence
56 the Pascal version could be useful in producing versions in other
57 structured languages.
59 Rather than use identifiers of minimal length (which therefore have
60 little mnemonic significance), the Pascal version uses meaningful
61 identifiers as follows [Note: A few changes have been made for C]:
64 BASIC C BASIC C BASIC C
66 A J S StickyBit
67 A1 AInverse J0 NoErrors T
68 B Radix [Failure] T0 Underflow
69 B1 BInverse J1 NoErrors T2 ThirtyTwo
70 B2 RadixD2 [SeriousDefect] T5 OneAndHalf
71 B9 BMinusU2 J2 NoErrors T7 TwentySeven
72 C [Defect] T8 TwoForty
73 C1 CInverse J3 NoErrors U OneUlp
74 D [Flaw] U0 UnderflowThreshold
75 D4 FourD K PageNo U1
76 E0 L Milestone U2
77 E1 M V
78 E2 Exp2 N V0
79 E3 N1 V8
80 E5 MinSqEr O Zero V9
81 E6 SqEr O1 One W
82 E7 MaxSqEr O2 Two X
83 E8 O3 Three X1
84 E9 O4 Four X8
85 F1 MinusOne O5 Five X9 Random1
86 F2 Half O8 Eight Y
87 F3 Third O9 Nine Y1
88 F6 P Precision Y2
89 F9 Q Y9 Random2
90 G1 GMult Q8 Z
91 G2 GDiv Q9 Z0 PseudoZero
92 G3 GAddSub R Z1
93 H R1 RMult Z2
94 H1 HInverse R2 RDiv Z9
95 I R3 RAddSub
96 IO NoTrials R4 RSqrt
97 I3 IEEE R9 Random9
99 SqRWrng
101 All the variables in BASIC are true variables and in consequence,
102 the program is more difficult to follow since the "constants" must
103 be determined (the glossary is very helpful). The Pascal version
104 uses Real constants, but checks are added to ensure that the values
105 are correctly converted by the compiler.
107 The major textual change to the Pascal version apart from the
108 identifiersis that named procedures are used, inserting parameters
109 wherehelpful. New procedures are also introduced. The
110 correspondence is as follows:
113 BASIC Pascal
114 lines
116 90- 140 Pause
117 170- 250 Instructions
118 380- 460 Heading
119 480- 670 Characteristics
120 690- 870 History
121 2940-2950 Random
122 3710-3740 NewD
123 4040-4080 DoesYequalX
124 4090-4110 PrintIfNPositive
125 4640-4850 TestPartialUnderflow
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.
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. */
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>
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. */
150 #define GTY(x)
152 #include "ansidecl.h"
153 #include "auto-host.h"
154 #include "hwint.h"
156 #undef EXTRA_MODES_FILE
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;
165 #define DEFTREECODE(SYM, STRING, TYPE, NARGS) SYM,
166 enum tree_code {
167 #include "tree.def"
168 LAST_AND_UNUSED_TREE_CODE
170 #undef DEFTREECODE
172 #define ENUM_BITFIELD(X) enum X
173 #define class klass
175 #include "real.h"
177 #undef class
180 /* We never produce signals from the library. Thus setjmp need do nothing. */
181 #undef setjmp
182 #define setjmp(x) (0)
184 static bool verbose = false;
185 static int verbose_index = 0;
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. */
192 class real_c_float
194 public:
195 static const enum machine_mode MODE = SFmode;
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];
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;
209 public:
210 real_c_float()
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)); }
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; }
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; }
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; }
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); }
247 const char * str () const;
248 const char * hex () const;
249 long integer () const;
250 int exp () const;
251 void ldexp (int);
254 void
255 real_c_float::from_long (long l)
257 REAL_VALUE_TYPE f;
259 real_from_integer (&f, MODE, l, l < 0 ? -1 : 0, 0);
260 real_to_target (image, &f, MODE);
263 void
264 real_c_float::from_str (const char *s)
266 REAL_VALUE_TYPE f;
267 const char *p = s;
269 if (*p == '-' || *p == '+')
270 p++;
271 if (strcasecmp(p, "inf") == 0)
273 real_inf (&f);
274 if (*s == '-')
275 real_arithmetic (&f, NEGATE_EXPR, &f, NULL);
277 else if (strcasecmp(p, "nan") == 0)
278 real_nan (&f, "", 1, MODE);
279 else
280 real_from_string (&f, s);
282 real_to_target (image, &f, MODE);
285 void
286 real_c_float::binop (int code, const real_c_float &b)
288 REAL_VALUE_TYPE ai, bi, ri;
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);
295 if (verbose)
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;
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);
307 switch (code)
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 ();
325 fprintf (stderr, "%6d: %s %c %s = %s\n", verbose_index++,
326 ab, symbol_for_code, bb, rb);
330 void
331 real_c_float::unop (int code)
333 REAL_VALUE_TYPE ai, ri;
335 real_from_target (&ai, image, MODE);
336 real_arithmetic (&ri, code, &ai, NULL);
337 real_to_target (image, &ri, MODE);
339 if (verbose)
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;
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);
350 switch (code)
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 ();
362 fprintf (stderr, "%6d: %s%s = %s\n", verbose_index++,
363 symbol_for_code, ab, rb);
367 bool
368 real_c_float::cmp (int code, const real_c_float &b) const
370 REAL_VALUE_TYPE ai, bi;
371 bool ret;
373 real_from_target (&ai, image, MODE);
374 real_from_target (&bi, b.image, MODE);
375 ret = real_compare (code, &ai, &bi);
377 if (verbose)
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;
384 real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
385 real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0);
387 switch (code)
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 ();
411 fprintf (stderr, "%6d: %s %s %s = %s\n", verbose_index++,
412 ab, symbol_for_code, bb, (ret ? "true" : "false"));
415 return ret;
418 const char *
419 real_c_float::str() const
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);
425 real_from_target (&f, image, MODE);
426 char *buf = new char[digits + 10];
427 real_to_decimal (buf, &f, digits+10, digits, 0);
429 return buf;
432 const char *
433 real_c_float::hex() const
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;
439 real_from_target (&f, image, MODE);
440 char *buf = new char[digits + 10];
441 real_to_hexadecimal (buf, &f, digits+10, digits, 0);
443 return buf;
446 long
447 real_c_float::integer() const
449 REAL_VALUE_TYPE f;
450 real_from_target (&f, image, MODE);
451 return real_to_integer (&f);
455 real_c_float::exp() const
457 REAL_VALUE_TYPE f;
458 real_from_target (&f, image, MODE);
459 return real_exponent (&f);
462 void
463 real_c_float::ldexp (int exp)
465 REAL_VALUE_TYPE ai;
467 real_from_target (&ai, image, MODE);
468 real_ldexp (&ai, &ai, exp);
469 real_to_target (image, &ai, MODE);
472 /* ====================================================================== */
473 /* An implementation of the abstract floating point class that uses native
474 arithmetic. Exists for reference and debugging. */
476 template<typename T>
477 class native_float
479 private:
480 // Force intermediate results back to memory.
481 volatile T image;
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);
489 public:
490 native_float()
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; }
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; }
506 const native_float& operator+= (const native_float &b)
508 image = verbose_binop(image, '+', b.image, image + b.image);
509 return *this;
511 const native_float& operator-= (const native_float &b)
513 image = verbose_binop(image, '-', b.image, image - b.image);
514 return *this;
516 const native_float& operator*= (const native_float &b)
518 image = verbose_binop(image, '*', b.image, image * b.image);
519 return *this;
521 const native_float& operator/= (const native_float &b)
523 image = verbose_binop(image, '/', b.image, image / b.image);
524 return *this;
527 native_float operator- () const
529 native_float r;
530 r.image = verbose_unop("-", image, -image);
531 return r;
533 native_float abs () const
535 native_float r;
536 r.image = verbose_unop("abs ", image, do_abs(image));
537 return r;
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); }
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);
561 template<typename T>
562 inline T
563 native_float<T>::from_str (const char *s)
565 return strtold (s, NULL);
568 template<>
569 inline float
570 native_float<float>::from_str (const char *s)
572 return strtof (s, NULL);
575 template<>
576 inline double
577 native_float<double>::from_str (const char *s)
579 return strtod (s, NULL);
582 template<typename T>
583 inline T
584 native_float<T>::do_abs (T image)
586 return fabsl (image);
589 template<>
590 inline float
591 native_float<float>::do_abs (float image)
593 return fabsf (image);
596 template<>
597 inline double
598 native_float<double>::do_abs (double image)
600 return fabs (image);
603 template<typename T>
605 native_float<T>::verbose_binop (T a, char symbol, T b, T r)
607 if (verbose)
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
620 return r;
623 template<typename T>
625 native_float<T>::verbose_unop (const char *symbol, T a, T r)
627 if (verbose)
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
638 return r;
641 template<typename T>
642 bool
643 native_float<T>::verbose_cmp (T a, const char *symbol, T b, bool r)
645 if (verbose)
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
658 return r;
661 template<typename T>
662 const char *
663 native_float<T>::str() const
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;
675 template<typename T>
676 const char *
677 native_float<T>::hex() const
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;
689 template<typename T>
691 native_float<T>::exp() const
693 int e;
694 frexp (image, &e);
695 return e;
698 template<typename T>
699 void
700 native_float<T>::ldexp (int exp)
702 image = ldexpl (image, exp);
705 template<>
706 void
707 native_float<float>::ldexp (int exp)
709 image = ldexpf (image, exp);
712 template<>
713 void
714 native_float<double>::ldexp (int exp)
716 image = ::ldexp (image, exp);
719 /* ====================================================================== */
720 /* Some libm routines that Paranoia expects to be available. */
722 template<typename FLOAT>
723 inline FLOAT
724 FABS (const FLOAT &f)
726 return f.abs();
729 template<typename FLOAT, typename RHS>
730 inline FLOAT
731 operator+ (const FLOAT &a, const RHS &b)
733 return FLOAT(a) += FLOAT(b);
736 template<typename FLOAT, typename RHS>
737 inline FLOAT
738 operator- (const FLOAT &a, const RHS &b)
740 return FLOAT(a) -= FLOAT(b);
743 template<typename FLOAT, typename RHS>
744 inline FLOAT
745 operator* (const FLOAT &a, const RHS &b)
747 return FLOAT(a) *= FLOAT(b);
750 template<typename FLOAT, typename RHS>
751 inline FLOAT
752 operator/ (const FLOAT &a, const RHS &b)
754 return FLOAT(a) /= FLOAT(b);
757 template<typename FLOAT>
758 FLOAT
759 FLOOR (const FLOAT &f)
761 /* ??? This is only correct when F is representable as an integer. */
762 long i = f.integer();
763 FLOAT r;
765 r = i;
766 if (i < 0 && f != r)
767 r = i - 1;
769 return r;
772 template<typename FLOAT>
773 FLOAT
774 SQRT (const FLOAT &f)
776 #if 0
777 FLOAT zero = long(0);
778 FLOAT two = 2;
779 FLOAT one = 1;
780 FLOAT diff, diff2;
781 FLOAT z, t;
783 if (f == zero)
784 return zero;
785 if (f < zero)
786 return zero / zero;
787 if (f == one)
788 return f;
790 z = f;
791 z.ldexp (-f.exp() / 2);
793 diff2 = FABS (z * z - f);
794 if (diff2 > zero)
795 while (1)
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;
805 return z;
806 #elif defined(NO_LONG_DOUBLE)
807 double d;
808 char buf[64];
810 d = strtod (f.hex(), NULL);
811 d = sqrt (d);
812 sprintf(buf, "%.35a", d);
814 return FLOAT(buf);
815 #else
816 long double ld;
817 char buf[64];
819 ld = strtold (f.hex(), NULL);
820 ld = sqrtl (ld);
821 sprintf(buf, "%.35La", ld);
823 return FLOAT(buf);
824 #endif
827 template<typename FLOAT>
828 FLOAT
829 LOG (FLOAT x)
831 #if 0
832 FLOAT zero = long(0);
833 FLOAT one = 1;
835 if (x <= zero)
836 return zero / zero;
837 if (x == one)
838 return zero;
840 int exp = x.exp() - 1;
841 x.ldexp(-exp);
843 FLOAT xm1 = x - one;
844 FLOAT y = xm1;
845 long n = 2;
847 FLOAT sum = xm1;
848 while (1)
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;
860 if (exp)
861 sum += FLOAT (exp) * FLOAT(".69314718055994530941");
863 return sum;
864 #elif defined (NO_LONG_DOUBLE)
865 double d;
866 char buf[64];
868 d = strtod (x.hex(), NULL);
869 d = log (d);
870 sprintf(buf, "%.35a", d);
872 return FLOAT(buf);
873 #else
874 long double ld;
875 char buf[64];
877 ld = strtold (x.hex(), NULL);
878 ld = logl (ld);
879 sprintf(buf, "%.35La", ld);
881 return FLOAT(buf);
882 #endif
885 template<typename FLOAT>
886 FLOAT
887 EXP (const FLOAT &x)
889 /* Cheat. */
890 #ifdef NO_LONG_DOUBLE
891 double d;
892 char buf[64];
894 d = strtod (x.hex(), NULL);
895 d = exp (d);
896 sprintf(buf, "%.35a", d);
898 return FLOAT(buf);
899 #else
900 long double ld;
901 char buf[64];
903 ld = strtold (x.hex(), NULL);
904 ld = expl (ld);
905 sprintf(buf, "%.35La", ld);
907 return FLOAT(buf);
908 #endif
911 template<typename FLOAT>
912 FLOAT
913 POW (const FLOAT &base, const FLOAT &exp)
915 /* Cheat. */
916 #ifdef NO_LONG_DOUBLE
917 double d1, d2;
918 char buf[64];
920 d1 = strtod (base.hex(), NULL);
921 d2 = strtod (exp.hex(), NULL);
922 d1 = pow (d1, d2);
923 sprintf(buf, "%.35a", d1);
925 return FLOAT(buf);
926 #else
927 long double ld1, ld2;
928 char buf[64];
930 ld1 = strtold (base.hex(), NULL);
931 ld2 = strtold (exp.hex(), NULL);
932 ld1 = powl (ld1, ld2);
933 sprintf(buf, "%.35La", ld1);
935 return FLOAT(buf);
936 #endif
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. */
943 int NoTrials = 20; /*Number of tests for commutativity. */
944 bool do_pause = false;
946 enum Guard { No, Yes };
947 enum Rounding { Other, Rounded, Chopped };
948 enum Class { Failure, Serious, Defect, Flaw };
950 template<typename FLOAT>
951 struct Paranoia
953 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
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;
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;
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 */
1007 int main ();
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 ();
1022 // Pretend we're bss.
1023 Paranoia() { memset(this, 0, sizeof (*this)); }
1026 template<typename FLOAT>
1028 Paranoia<FLOAT>::main()
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");
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");
1060 Z = -Zero;
1061 if (Z != Zero)
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 ();
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");
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");
1082 TstCond (Failure, Half + MinusOne + Half == Zero, "1/2 + (-1) + 1/2 != 0");
1084 /*=============================================*/
1085 Milestone = 10;
1086 /*=============================================*/
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");
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");
1104 if (ErrCnt[Failure] == 0)
1106 printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
1107 printf ("\n");
1109 printf ("Searching for Radix and Precision.\n");
1110 W = One;
1113 W = W + W;
1114 Y = W + One;
1115 Z = Y - W;
1116 Y = Z - One;
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;
1124 Radix = W + Y;
1125 Y = Y + Y;
1126 Radix = Radix - W;
1128 while (Radix == Zero);
1129 if (Radix < Two)
1130 Radix = One;
1131 printf ("Radix = %s .\n", Radix.str());
1132 if (Radix != One)
1134 W = One;
1137 Precision = Precision + One;
1138 W = W * Radix;
1139 Y = W + One;
1141 while ((Y - W) == One);
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 ");
1150 /*save old values */
1151 E0 = Radix;
1152 E1 = U1;
1153 E9 = U2;
1154 E3 = Precision;
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;
1164 /*... now X = (unknown no.) ulps of 1+... */
1167 U2 = X;
1168 Y = Half * U2 + ThirtyTwo * U2 * U2;
1169 Y = One + Y;
1170 X = Y - One;
1172 while (!((U2 <= X) || (X <= Zero)));
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;
1183 /*... now X == (unknown no.) ulps of 1 -... */
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;
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;
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)
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;
1242 if ((Precision != FLOOR (Precision)) || (Radix == One))
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");
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);
1263 Z2 = X;
1264 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
1266 while (!((Z2 <= X) || (X <= Zero)));
1267 X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four);
1270 Z1 = Z;
1271 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
1272 + One / Two)) + One / Two;
1274 while (!((Z1 <= Z) || (Z <= Zero)));
1279 Y1 = Y;
1281 (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half)) +
1282 Half;
1284 while (!((Y1 <= Y) || (Y <= Zero)));
1285 X1 = X;
1286 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
1288 while (!((X1 <= X) || (X <= Zero)));
1289 if ((X1 != Y1) || (X1 != Z1))
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");
1299 else
1301 if ((Z1 != U1) || (Z2 != U2))
1303 if ((Z1 >= U1) || (Z2 >= U2))
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());
1310 else
1312 if ((Z1 <= Zero) || (Z2 <= Zero))
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");
1319 if (Z1 != Z2 || Z1 > Zero)
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());
1333 printf
1334 ("That feature is not tested further by this program.\n");
1338 Pause ();
1339 /*=============================================*/
1340 Milestone = 35;
1341 /*=============================================*/
1342 if (Radix >= Two)
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.");
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
1373 GMult = No;
1374 TstCond (Serious, false, "* lacks a Guard Digit, so 1*X != X");
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
1405 GDiv = No;
1406 TstCond (Serious, false, "Division lacks a Guard Digit, so X/1 != X");
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
1431 GAddSub = No;
1432 TstCond (Serious, false,
1433 "- lacks Guard Digit, so cancellation is obscured");
1435 if (F9 != One && F9 - One >= Zero)
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");
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;
1458 AInvrse = Radix;
1461 X = AInvrse;
1462 AInvrse = AInvrse / A1;
1464 while (!(FLOOR (AInvrse) != AInvrse));
1465 Done = (X == One) || (A1 > Three);
1466 if (!Done)
1467 A1 = Nine + One;
1469 while (!(Done));
1470 if (X == One)
1471 A1 = Radix;
1472 AInvrse = One / A1;
1473 X = A1;
1474 Y = AInvrse;
1475 Done = false;
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;
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))
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))
1516 RMult = Rounded;
1517 printf ("Multiplication appears to round correctly.\n");
1519 else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
1520 && (T < Zero) && (StickyBit + U2 == Zero) && (Y1 < Half))
1522 RMult = Chopped;
1523 printf ("Multiplication appears to chop.\n");
1525 else
1526 printf ("* is neither chopped nor correctly rounded.\n");
1527 if ((RMult == Rounded) && (GMult == No))
1528 notify ("Multiplication");
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)))
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))
1565 RDiv = Rounded;
1566 printf ("Division appears to round correctly.\n");
1567 if (GDiv == No)
1568 notify ("Division");
1570 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
1571 && (Y2 < Zero) && (Y1 - Half < F9 - Half))
1573 RDiv = Chopped;
1574 printf ("Division appears to chop.\n");
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))
1595 RAddSub = Chopped;
1596 printf ("Add/Subtract appears to be chopped.\n");
1598 if (GAddSub == Yes)
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))
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))
1616 RAddSub = Rounded;
1617 printf ("Addition/Subtraction appears to round correctly.\n");
1618 if (GAddSub == No)
1619 notify ("Add/Subtract");
1621 else
1622 printf ("Addition/Subtraction neither rounds nor chops.\n");
1624 else
1625 printf ("Addition/Subtraction neither rounds nor chops.\n");
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)
1637 S = Zero;
1638 BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n");
1640 StickyBit = Zero;
1641 if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
1642 && (RMult == Rounded) && (RDiv == Rounded)
1643 && (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2))
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))
1652 Z = T + Y;
1653 Y = Z - X;
1654 if ((Z - T >= U2) && (Y - T == Zero))
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))
1662 Z = (Half - U1) * U1;
1663 T = F9 - Z;
1664 Q = F9 - Y;
1665 if ((T - F9 == Zero) && (F9 - U1 - Q == Zero))
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)
1674 if (Radix != Two)
1676 X = Two + U2;
1677 Y = X / Two;
1678 if ((Y - One == Zero))
1679 StickyBit = S;
1681 else
1682 StickyBit = S;
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;
1708 X = Random ();
1709 Y = Random ();
1710 Z9 = Y * X;
1711 Z = X * Y;
1712 Z9 = Z - Z9;
1713 I = I + 1;
1715 while (!((I > NoTrials) || (Z9 != Zero)));
1716 if (I == NoTrials)
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);
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))
1758 X = Y;
1759 Y = Radix * Y;
1761 while (!((Y - X >= NoTrials)));
1762 OneUlp = X * U2;
1763 I = 1;
1764 while (I <= NoTrials)
1766 X = X + One;
1767 SqXMinX (Defect);
1768 if (J > Zero)
1769 break;
1770 I = I + 1;
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))
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
1789 Q = FLOOR (Q + Half);
1790 if (!(I > 0 || Radix == Q * Q))
1791 Monot = true;
1792 else if (I > 0)
1794 if (I > 1)
1795 Monot = true;
1796 else
1798 Y = Y * BInvrse;
1799 X = Y - U1;
1800 Z = Y + U1;
1803 else
1805 Y = Q;
1806 X = Y - U2;
1807 Z = Y + U2;
1811 if (Monot)
1812 printf ("sqrt has passed a test for Monotonicity.\n");
1813 else
1815 BadCond (Defect, "");
1816 printf ("sqrt(X) is non-monotonic for X near %s .\n", Y.str());
1818 /*=============================================*/
1819 Milestone = 110;
1820 /*=============================================*/
1821 printf ("Seeking Underflow thresholds UfThold and E0.\n");
1822 D = U1;
1823 if (Precision != FLOOR (Precision))
1825 D = BInvrse;
1826 X = Precision;
1829 D = D * BInvrse;
1830 X = X - One;
1832 while (X > Zero);
1834 Y = One;
1835 Z = D;
1836 /* ... D is power of 1/Radix < 1. */
1839 C = Y;
1840 Y = Z;
1841 Z = Y * Y;
1843 while ((Y > Z) && (Z + Z > Z));
1844 Y = C;
1845 Z = Y * D;
1848 C = Y;
1849 Y = Z;
1850 Z = Y * D;
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 */
1865 Y = E0;
1866 E0 = Z;
1867 Z = E0 * H;
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)
1878 E9 = Radix * U2;
1879 S = One + E9;
1880 D = C * S;
1881 if (D <= C)
1883 BadCond (Failure,
1884 "multiplication gets too many last digits wrong.\n");
1885 Underflow = E0;
1886 Y1 = Zero;
1887 PseudoZero = Z;
1888 Pause ();
1891 else
1893 Underflow = D;
1894 PseudoZero = Underflow * H;
1895 UfThold = Zero;
1898 Y1 = Underflow;
1899 Underflow = PseudoZero;
1900 if (E1 + E1 <= E1)
1902 Y2 = Underflow * HInvrse;
1903 E1 = FABS (Y1 - Y2);
1904 Q = Y1;
1905 if ((UfThold == Zero) && (Y1 != Y2))
1906 UfThold = Y1;
1908 PseudoZero = PseudoZero * H;
1910 while ((Underflow > PseudoZero)
1911 && (PseudoZero + PseudoZero > PseudoZero));
1913 /* Comment line 4530 .. 4560 */
1914 if (PseudoZero != Zero)
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)
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)
1929 printf ("But -PseudoZero, which should be\n");
1930 printf ("positive, isn't; it prints out as %s .\n", X.str());
1933 else
1935 BadCond (Flaw, "Underflow can stick at an allegedly positive\n");
1936 printf ("value PseudoZero that prints out as %s .\n",
1937 PseudoZero.str());
1939 TstPtUf ();
1941 /*=============================================*/
1942 Milestone = 120;
1943 /*=============================================*/
1944 if (CInvrse * Y > CInvrse * Y1)
1946 S = H * S;
1947 E0 = Underflow;
1949 if (!((E1 == Zero) || (E1 == E0)))
1951 BadCond (Defect, "");
1952 if (E1 < E0)
1954 printf ("Products underflow at a higher");
1955 printf (" threshold than differences.\n");
1956 if (PseudoZero == Zero)
1957 E0 = E1;
1959 else
1961 printf ("Difference underflows at a higher");
1962 printf (" threshold than products.\n");
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)
1979 case 1:
1980 UfThold = Underflow;
1981 if ((CInvrse * Q) != ((CInvrse * Y) * S))
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");
1990 Pause ();
1991 break;
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;
2002 case 3:
2003 X = X;
2004 break;
2006 case 4:
2007 if ((Q == UfThold) && (E1 == E0) && (FABS (UfThold - E1 / E9) <= E1))
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);
2019 if (UfNGrad)
2021 printf ("\n");
2022 if (setjmp (ovfl_buf))
2024 printf ("Underflow / UfThold failed!\n");
2025 R = H + H;
2027 else
2028 R = SQRT (Underflow / UfThold);
2029 if (R <= H)
2031 Z = R * UfThold;
2032 X = Z * (One + R * H * (One + H));
2034 else
2036 Z = UfThold;
2037 X = Z * (One + H * H * (One + H));
2039 if (!((X == Z) || (X - Z != Zero)))
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());
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)
2065 if (Y > E0)
2067 BadCond (Defect, "");
2068 I = 5;
2070 else
2072 BadCond (Serious, "");
2073 I = 4;
2075 printf ("Range is too narrow; U1^%d Underflows.\n", I);
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))
2089 BadCond (Serious, "trap on underflow.\n");
2091 else
2093 V9 = POW (HInvrse, Y2);
2094 printf (" %s .\n", V9.str());
2095 if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold)))
2097 BadCond (Serious, "this is not between 0 and underflow\n");
2098 printf (" threshold = %s .\n", UfThold.str());
2100 else if (!(V9 > UfThold * (One + E9)))
2101 printf ("This computed value is O.K.\n");
2102 else
2104 BadCond (Defect, "this is not between 0 and underflow\n");
2105 printf (" threshold = %s .\n", UfThold.str());
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))
2118 I = 0;
2119 V9 = Y;
2120 goto overflow;
2124 V = Y;
2125 Y = V9;
2126 V9 = HInvrse * Y;
2128 while (V9 < Y);
2129 I = 1;
2130 overflow:
2131 Z = V9;
2132 printf ("Can `Z = -Y' overflow?\n");
2133 printf ("Trying it on Y = %s .\n", Y.str());
2134 V9 = -Y;
2135 V0 = V9;
2136 if (V - Y == V + V0)
2137 printf ("Seems O.K.\n");
2138 else
2140 printf ("finds a ");
2141 BadCond (Flaw, "-(-Y) differs from Y.\n");
2143 if (Z != Y)
2145 BadCond (Serious, "");
2146 printf ("overflow past %s\n\tshrinks to %s .\n", Y.str(), Z.str());
2148 if (I)
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;
2159 else
2161 V = Y * (HInvrse * U2 - HInvrse);
2162 V = V + ((One - HInvrse) * U2) * Y;
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))
2181 BadCond (Failure, "Comparisons involving ");
2182 printf ("+-%s, +-%s\nand +-%s are confused by Overflow.",
2183 V.str(), V0.str(), UfThold.str());
2185 /*=============================================*/
2186 Milestone = 175;
2187 /*=============================================*/
2188 printf ("\n");
2189 for (Indx = 1; Indx <= 3; ++Indx)
2191 switch (Indx)
2193 case 1:
2194 Z = UfThold;
2195 break;
2196 case 2:
2197 Z = E0;
2198 break;
2199 case 3:
2200 Z = PseudoZero;
2201 break;
2203 if (Z != Zero)
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());
2219 /*=============================================*/
2220 Milestone = 180;
2221 /*=============================================*/
2222 for (Indx = 1; Indx <= 2; ++Indx)
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)))
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());
2242 /*=============================================*/
2243 Milestone = 190;
2244 /*=============================================*/
2245 Pause ();
2246 X = UfThold * V;
2247 Y = Radix * Radix;
2248 if (X * Y < One || X > Y)
2250 if (X * Y < U1 || X > Y / U1)
2251 BadCond (Defect, "Badly");
2252 else
2253 BadCond (Flaw, "");
2255 printf (" unbalanced range; UfThold * V = %s\n\t%s\n",
2256 X.str(), "is too far from 1.\n");
2258 /*=============================================*/
2259 Milestone = 200;
2260 /*=============================================*/
2261 for (Indx = 1; Indx <= 5; ++Indx)
2263 X = F9;
2264 switch (Indx)
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;
2278 Y = X;
2279 if (setjmp (ovfl_buf))
2280 printf (" X / X traps when X = %s\n", X.str());
2281 else
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());
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");
2312 static const char *msg[] = {
2313 "FAILUREs encountered =",
2314 "SERIOUS DEFECTs discovered =",
2315 "DEFECTs discovered =",
2316 "FLAWs discovered ="
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]);
2323 printf ("\n");
2324 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] + ErrCnt[Flaw]) > 0)
2326 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] == 0)
2327 && (ErrCnt[Flaw] > 0))
2329 printf ("The arithmetic diagnosed seems ");
2330 printf ("Satisfactory though flawed.\n");
2332 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0) && (ErrCnt[Defect] > 0))
2334 printf ("The arithmetic diagnosed may be Acceptable\n");
2335 printf ("despite inconvenient Defects.\n");
2337 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0)
2339 printf ("The arithmetic diagnosed has ");
2340 printf ("unacceptable Serious Defects.\n");
2342 if (ErrCnt[Failure] > 0)
2344 printf ("Potentially fatal FAILURE may have spoiled this");
2345 printf (" program's subsequent diagnoses.\n");
2348 else
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
2356 if (StickyBit >= One &&
2357 (Radix - Two) * (Radix - Nine - One) == Zero)
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
2371 printf (",\nexcept for possibly Double Rounding");
2372 printf (" during Gradual Underflow.\n");
2375 printf ("The arithmetic diagnosed appears to be Excellent!\n");
2378 printf ("END OF TEST.\n");
2379 return 0;
2382 template<typename FLOAT>
2383 FLOAT
2384 Paranoia<FLOAT>::Sign (FLOAT X)
2386 return X >= FLOAT (long (0)) ? 1 : -1;
2389 template<typename FLOAT>
2390 void
2391 Paranoia<FLOAT>::Pause ()
2393 if (do_pause)
2395 fputs ("Press return...", stdout);
2396 fflush (stdout);
2397 getchar();
2399 printf ("\nDiagnosis resumes after milestone Number %d", Milestone);
2400 printf (" Page: %d\n\n", PageNo);
2401 ++Milestone;
2402 ++PageNo;
2405 template<typename FLOAT>
2406 void
2407 Paranoia<FLOAT>::TstCond (int K, int Valid, const char *T)
2409 if (!Valid)
2411 BadCond (K, T);
2412 printf (".\n");
2416 template<typename FLOAT>
2417 void
2418 Paranoia<FLOAT>::BadCond (int K, const char *T)
2420 static const char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
2422 ErrCnt[K] = ErrCnt[K] + 1;
2423 printf ("%s: %s", msg[K], T);
2426 /* Random computes
2427 X = (Random1 + Random9)^5
2428 Random1 = X - FLOOR(X) + 0.000005 * X;
2429 and returns the new value of Random1. */
2431 template<typename FLOAT>
2432 FLOAT
2433 Paranoia<FLOAT>::Random ()
2435 FLOAT X, Y;
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);
2446 template<typename FLOAT>
2447 void
2448 Paranoia<FLOAT>::SqXMinX (int ErrKind)
2450 FLOAT XA, XB;
2452 XB = X * BInvrse;
2453 XA = X - XB;
2454 SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp;
2455 if (SqEr != Zero)
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");
2469 template<typename FLOAT>
2470 void
2471 Paranoia<FLOAT>::NewD ()
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)
2479 Z = -Z;
2480 Z1 = -Z1;
2482 D = Radix * D;
2485 template<typename FLOAT>
2486 void
2487 Paranoia<FLOAT>::SR3750 ()
2489 if (!((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2)))
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;
2505 template<typename FLOAT>
2506 void
2507 Paranoia<FLOAT>::IsYeqX ()
2509 if (Y != X)
2511 if (N <= 0)
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());
2522 N = N + 1; /* ... count discrepancies. */
2526 template<typename FLOAT>
2527 void
2528 Paranoia<FLOAT>::PrintIfNPositive ()
2530 if (N > 0)
2531 printf ("Similar discrepancies have occurred %d times.\n", N);
2534 template<typename FLOAT>
2535 void
2536 Paranoia<FLOAT>::TstPtUf ()
2538 N = 0;
2539 if (Z != Zero)
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)
2549 printf ("This is O.K., provided Over/Underflow");
2550 printf (" has NOT just been signaled.\n");
2552 else
2554 if ((Q9 < One) || (Q9 > Two))
2556 very_serious:
2557 N = 1;
2558 ErrCnt[Serious] = ErrCnt[Serious] + 1;
2559 printf ("This is a VERY SERIOUS DEFECT!\n");
2561 else
2563 N = 1;
2564 ErrCnt[Defect] = ErrCnt[Defect] + 1;
2565 printf ("This is a DEFECT!\n");
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))
2575 if (N > 0)
2576 Pause ();
2578 else
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)
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());
2596 Pause ();
2601 template<typename FLOAT>
2602 void
2603 Paranoia<FLOAT>::notify (const char *s)
2605 printf ("%s test appears to be inconsistent...\n", s);
2606 printf (" PLEASE NOTIFY KARPINKSI!\n");
2609 /* ====================================================================== */
2611 int main(int ac, char **av)
2613 setbuf(stdout, NULL);
2614 setbuf(stderr, NULL);
2616 while (1)
2617 switch (getopt (ac, av, "pvg:fdl"))
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':
2629 static const struct {
2630 const char *name;
2631 const struct real_format *fmt;
2632 } fmts[] = {
2633 #define F(x) { #x, &x##_format }
2634 F(ieee_single),
2635 F(ieee_double),
2636 F(ieee_extended_motorola),
2637 F(ieee_extended_intel_96),
2638 F(ieee_extended_intel_128),
2639 F(ibm_extended),
2640 F(ieee_quad),
2641 F(vax_f),
2642 F(vax_d),
2643 F(vax_g),
2644 F(i370_single),
2645 F(i370_double),
2646 F(c4x_single),
2647 F(c4x_extended),
2648 F(real_internal),
2649 #undef F
2652 int i, n = sizeof (fmts)/sizeof(*fmts);
2654 for (i = 0; i < n; ++i)
2655 if (strcmp (fmts[i].name, optarg) == 0)
2656 break;
2658 if (i == n)
2660 printf ("Unknown implementation \"%s\"; "
2661 "available implementations:\n", optarg);
2662 for (i = 0; i < n; ++i)
2663 printf ("\t%s\n", fmts[i].name);
2664 return 1;
2667 // We cheat and use the same mode all the time, but vary
2668 // the format used for that mode.
2669 real_format_for_mode[int(real_c_float::MODE) - int(QFmode)]
2670 = fmts[i].fmt;
2672 Paranoia<real_c_float>().main();
2673 break;
2676 case 'f':
2677 Paranoia < native_float<float> >().main();
2678 break;
2679 case 'd':
2680 Paranoia < native_float<double> >().main();
2681 break;
2682 case 'l':
2683 #ifndef NO_LONG_DOUBLE
2684 Paranoia < native_float<long double> >().main();
2685 #endif
2686 break;
2688 case '?':
2689 puts ("-p\tpause between pages");
2690 puts ("-g<FMT>\treal.c implementation FMT");
2691 puts ("-f\tnative float");
2692 puts ("-d\tnative double");
2693 puts ("-l\tnative long double");
2694 return 0;
2698 /* GCC stuff referenced by real.o. */
2700 extern "C" void
2701 fancy_abort ()
2703 abort ();
2706 int target_flags = 0;
2708 extern "C" int
2709 floor_log2_wide (unsigned HOST_WIDE_INT x)
2711 int log = -1;
2712 while (x != 0)
2713 log++,
2714 x >>= 1;
2715 return log;