* vi.po: Update.
[official-gcc.git] / contrib / paranoia.cc
blob8e8500e23898dbcb290f60864a3729eba46bad4d
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 class klass
174 #include "real.h"
176 #undef class
179 /* We never produce signals from the library. Thus setjmp need do nothing. */
180 #undef setjmp
181 #define setjmp(x) (0)
183 static bool verbose = false;
184 static int verbose_index = 0;
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. */
191 class real_c_float
193 public:
194 static const enum machine_mode MODE = SFmode;
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];
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;
208 public:
209 real_c_float()
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)); }
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; }
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; }
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; }
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); }
246 const char * str () const;
247 const char * hex () const;
248 long integer () const;
249 int exp () const;
250 void ldexp (int);
253 void
254 real_c_float::from_long (long l)
256 REAL_VALUE_TYPE f;
258 real_from_integer (&f, MODE, l, l < 0 ? -1 : 0, 0);
259 real_to_target (image, &f, MODE);
262 void
263 real_c_float::from_str (const char *s)
265 REAL_VALUE_TYPE f;
266 const char *p = s;
268 if (*p == '-' || *p == '+')
269 p++;
270 if (strcasecmp(p, "inf") == 0)
272 real_inf (&f);
273 if (*s == '-')
274 real_arithmetic (&f, NEGATE_EXPR, &f, NULL);
276 else if (strcasecmp(p, "nan") == 0)
277 real_nan (&f, "", 1, MODE);
278 else
279 real_from_string (&f, s);
281 real_to_target (image, &f, MODE);
284 void
285 real_c_float::binop (int code, const real_c_float &b)
287 REAL_VALUE_TYPE ai, bi, ri;
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);
294 if (verbose)
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;
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);
306 switch (code)
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 ();
324 fprintf (stderr, "%6d: %s %c %s = %s\n", verbose_index++,
325 ab, symbol_for_code, bb, rb);
329 void
330 real_c_float::unop (int code)
332 REAL_VALUE_TYPE ai, ri;
334 real_from_target (&ai, image, MODE);
335 real_arithmetic (&ri, code, &ai, NULL);
336 real_to_target (image, &ri, MODE);
338 if (verbose)
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;
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);
349 switch (code)
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 ();
361 fprintf (stderr, "%6d: %s%s = %s\n", verbose_index++,
362 symbol_for_code, ab, rb);
366 bool
367 real_c_float::cmp (int code, const real_c_float &b) const
369 REAL_VALUE_TYPE ai, bi;
370 bool ret;
372 real_from_target (&ai, image, MODE);
373 real_from_target (&bi, b.image, MODE);
374 ret = real_compare (code, &ai, &bi);
376 if (verbose)
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;
383 real_to_hexadecimal (ab, &ai, sizeof(ab), digits, 0);
384 real_to_hexadecimal (bb, &bi, sizeof(bb), digits, 0);
386 switch (code)
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 ();
410 fprintf (stderr, "%6d: %s %s %s = %s\n", verbose_index++,
411 ab, symbol_for_code, bb, (ret ? "true" : "false"));
414 return ret;
417 const char *
418 real_c_float::str() const
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);
424 real_from_target (&f, image, MODE);
425 char *buf = new char[digits + 10];
426 real_to_decimal (buf, &f, digits+10, digits, 0);
428 return buf;
431 const char *
432 real_c_float::hex() const
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;
438 real_from_target (&f, image, MODE);
439 char *buf = new char[digits + 10];
440 real_to_hexadecimal (buf, &f, digits+10, digits, 0);
442 return buf;
445 long
446 real_c_float::integer() const
448 REAL_VALUE_TYPE f;
449 real_from_target (&f, image, MODE);
450 return real_to_integer (&f);
454 real_c_float::exp() const
456 REAL_VALUE_TYPE f;
457 real_from_target (&f, image, MODE);
458 return real_exponent (&f);
461 void
462 real_c_float::ldexp (int exp)
464 REAL_VALUE_TYPE ai;
466 real_from_target (&ai, image, MODE);
467 real_ldexp (&ai, &ai, exp);
468 real_to_target (image, &ai, MODE);
471 /* ====================================================================== */
472 /* An implementation of the abstract floating point class that uses native
473 arithmetic. Exists for reference and debugging. */
475 template<typename T>
476 class native_float
478 private:
479 // Force intermediate results back to memory.
480 volatile T image;
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);
488 public:
489 native_float()
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; }
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; }
505 const native_float& operator+= (const native_float &b)
507 image = verbose_binop(image, '+', b.image, image + b.image);
508 return *this;
510 const native_float& operator-= (const native_float &b)
512 image = verbose_binop(image, '-', b.image, image - b.image);
513 return *this;
515 const native_float& operator*= (const native_float &b)
517 image = verbose_binop(image, '*', b.image, image * b.image);
518 return *this;
520 const native_float& operator/= (const native_float &b)
522 image = verbose_binop(image, '/', b.image, image / b.image);
523 return *this;
526 native_float operator- () const
528 native_float r;
529 r.image = verbose_unop("-", image, -image);
530 return r;
532 native_float abs () const
534 native_float r;
535 r.image = verbose_unop("abs ", image, do_abs(image));
536 return r;
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); }
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);
560 template<typename T>
561 inline T
562 native_float<T>::from_str (const char *s)
564 return strtold (s, NULL);
567 template<>
568 inline float
569 native_float<float>::from_str (const char *s)
571 return strtof (s, NULL);
574 template<>
575 inline double
576 native_float<double>::from_str (const char *s)
578 return strtod (s, NULL);
581 template<typename T>
582 inline T
583 native_float<T>::do_abs (T image)
585 return fabsl (image);
588 template<>
589 inline float
590 native_float<float>::do_abs (float image)
592 return fabsf (image);
595 template<>
596 inline double
597 native_float<double>::do_abs (double image)
599 return fabs (image);
602 template<typename T>
604 native_float<T>::verbose_binop (T a, char symbol, T b, T r)
606 if (verbose)
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
619 return r;
622 template<typename T>
624 native_float<T>::verbose_unop (const char *symbol, T a, T r)
626 if (verbose)
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
637 return r;
640 template<typename T>
641 bool
642 native_float<T>::verbose_cmp (T a, const char *symbol, T b, bool r)
644 if (verbose)
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
657 return r;
660 template<typename T>
661 const char *
662 native_float<T>::str() const
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;
674 template<typename T>
675 const char *
676 native_float<T>::hex() const
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;
688 template<typename T>
690 native_float<T>::exp() const
692 int e;
693 frexp (image, &e);
694 return e;
697 template<typename T>
698 void
699 native_float<T>::ldexp (int exp)
701 image = ldexpl (image, exp);
704 template<>
705 void
706 native_float<float>::ldexp (int exp)
708 image = ldexpf (image, exp);
711 template<>
712 void
713 native_float<double>::ldexp (int exp)
715 image = ::ldexp (image, exp);
718 /* ====================================================================== */
719 /* Some libm routines that Paranoia expects to be available. */
721 template<typename FLOAT>
722 inline FLOAT
723 FABS (const FLOAT &f)
725 return f.abs();
728 template<typename FLOAT, typename RHS>
729 inline FLOAT
730 operator+ (const FLOAT &a, const RHS &b)
732 return FLOAT(a) += FLOAT(b);
735 template<typename FLOAT, typename RHS>
736 inline FLOAT
737 operator- (const FLOAT &a, const RHS &b)
739 return FLOAT(a) -= FLOAT(b);
742 template<typename FLOAT, typename RHS>
743 inline FLOAT
744 operator* (const FLOAT &a, const RHS &b)
746 return FLOAT(a) *= FLOAT(b);
749 template<typename FLOAT, typename RHS>
750 inline FLOAT
751 operator/ (const FLOAT &a, const RHS &b)
753 return FLOAT(a) /= FLOAT(b);
756 template<typename FLOAT>
757 FLOAT
758 FLOOR (const FLOAT &f)
760 /* ??? This is only correct when F is representable as an integer. */
761 long i = f.integer();
762 FLOAT r;
764 r = i;
765 if (i < 0 && f != r)
766 r = i - 1;
768 return r;
771 template<typename FLOAT>
772 FLOAT
773 SQRT (const FLOAT &f)
775 #if 0
776 FLOAT zero = long(0);
777 FLOAT two = 2;
778 FLOAT one = 1;
779 FLOAT diff, diff2;
780 FLOAT z, t;
782 if (f == zero)
783 return zero;
784 if (f < zero)
785 return zero / zero;
786 if (f == one)
787 return f;
789 z = f;
790 z.ldexp (-f.exp() / 2);
792 diff2 = FABS (z * z - f);
793 if (diff2 > zero)
794 while (1)
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;
804 return z;
805 #elif defined(NO_LONG_DOUBLE)
806 double d;
807 char buf[64];
809 d = strtod (f.hex(), NULL);
810 d = sqrt (d);
811 sprintf(buf, "%.35a", d);
813 return FLOAT(buf);
814 #else
815 long double ld;
816 char buf[64];
818 ld = strtold (f.hex(), NULL);
819 ld = sqrtl (ld);
820 sprintf(buf, "%.35La", ld);
822 return FLOAT(buf);
823 #endif
826 template<typename FLOAT>
827 FLOAT
828 LOG (FLOAT x)
830 #if 0
831 FLOAT zero = long(0);
832 FLOAT one = 1;
834 if (x <= zero)
835 return zero / zero;
836 if (x == one)
837 return zero;
839 int exp = x.exp() - 1;
840 x.ldexp(-exp);
842 FLOAT xm1 = x - one;
843 FLOAT y = xm1;
844 long n = 2;
846 FLOAT sum = xm1;
847 while (1)
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;
859 if (exp)
860 sum += FLOAT (exp) * FLOAT(".69314718055994530941");
862 return sum;
863 #elif defined (NO_LONG_DOUBLE)
864 double d;
865 char buf[64];
867 d = strtod (x.hex(), NULL);
868 d = log (d);
869 sprintf(buf, "%.35a", d);
871 return FLOAT(buf);
872 #else
873 long double ld;
874 char buf[64];
876 ld = strtold (x.hex(), NULL);
877 ld = logl (ld);
878 sprintf(buf, "%.35La", ld);
880 return FLOAT(buf);
881 #endif
884 template<typename FLOAT>
885 FLOAT
886 EXP (const FLOAT &x)
888 /* Cheat. */
889 #ifdef NO_LONG_DOUBLE
890 double d;
891 char buf[64];
893 d = strtod (x.hex(), NULL);
894 d = exp (d);
895 sprintf(buf, "%.35a", d);
897 return FLOAT(buf);
898 #else
899 long double ld;
900 char buf[64];
902 ld = strtold (x.hex(), NULL);
903 ld = expl (ld);
904 sprintf(buf, "%.35La", ld);
906 return FLOAT(buf);
907 #endif
910 template<typename FLOAT>
911 FLOAT
912 POW (const FLOAT &base, const FLOAT &exp)
914 /* Cheat. */
915 #ifdef NO_LONG_DOUBLE
916 double d1, d2;
917 char buf[64];
919 d1 = strtod (base.hex(), NULL);
920 d2 = strtod (exp.hex(), NULL);
921 d1 = pow (d1, d2);
922 sprintf(buf, "%.35a", d1);
924 return FLOAT(buf);
925 #else
926 long double ld1, ld2;
927 char buf[64];
929 ld1 = strtold (base.hex(), NULL);
930 ld2 = strtold (exp.hex(), NULL);
931 ld1 = powl (ld1, ld2);
932 sprintf(buf, "%.35La", ld1);
934 return FLOAT(buf);
935 #endif
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. */
942 int NoTrials = 20; /*Number of tests for commutativity. */
943 bool do_pause = false;
945 enum Guard { No, Yes };
946 enum Rounding { Other, Rounded, Chopped };
947 enum Class { Failure, Serious, Defect, Flaw };
949 template<typename FLOAT>
950 struct Paranoia
952 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
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;
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;
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 */
1006 int main ();
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 ();
1021 // Pretend we're bss.
1022 Paranoia() { memset(this, 0, sizeof (*this)); }
1025 template<typename FLOAT>
1027 Paranoia<FLOAT>::main()
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");
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");
1059 Z = -Zero;
1060 if (Z != Zero)
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 ();
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");
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");
1081 TstCond (Failure, Half + MinusOne + Half == Zero, "1/2 + (-1) + 1/2 != 0");
1083 /*=============================================*/
1084 Milestone = 10;
1085 /*=============================================*/
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");
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");
1103 if (ErrCnt[Failure] == 0)
1105 printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
1106 printf ("\n");
1108 printf ("Searching for Radix and Precision.\n");
1109 W = One;
1112 W = W + W;
1113 Y = W + One;
1114 Z = Y - W;
1115 Y = Z - One;
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;
1123 Radix = W + Y;
1124 Y = Y + Y;
1125 Radix = Radix - W;
1127 while (Radix == Zero);
1128 if (Radix < Two)
1129 Radix = One;
1130 printf ("Radix = %s .\n", Radix.str());
1131 if (Radix != One)
1133 W = One;
1136 Precision = Precision + One;
1137 W = W * Radix;
1138 Y = W + One;
1140 while ((Y - W) == One);
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 ");
1149 /*save old values */
1150 E0 = Radix;
1151 E1 = U1;
1152 E9 = U2;
1153 E3 = Precision;
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;
1163 /*... now X = (unknown no.) ulps of 1+... */
1166 U2 = X;
1167 Y = Half * U2 + ThirtyTwo * U2 * U2;
1168 Y = One + Y;
1169 X = Y - One;
1171 while (!((U2 <= X) || (X <= Zero)));
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;
1182 /*... now X == (unknown no.) ulps of 1 -... */
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;
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;
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)
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;
1241 if ((Precision != FLOOR (Precision)) || (Radix == One))
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");
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);
1262 Z2 = X;
1263 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
1265 while (!((Z2 <= X) || (X <= Zero)));
1266 X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four);
1269 Z1 = Z;
1270 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
1271 + One / Two)) + One / Two;
1273 while (!((Z1 <= Z) || (Z <= Zero)));
1278 Y1 = Y;
1280 (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half)) +
1281 Half;
1283 while (!((Y1 <= Y) || (Y <= Zero)));
1284 X1 = X;
1285 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
1287 while (!((X1 <= X) || (X <= Zero)));
1288 if ((X1 != Y1) || (X1 != Z1))
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");
1298 else
1300 if ((Z1 != U1) || (Z2 != U2))
1302 if ((Z1 >= U1) || (Z2 >= U2))
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());
1309 else
1311 if ((Z1 <= Zero) || (Z2 <= Zero))
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");
1318 if (Z1 != Z2 || Z1 > Zero)
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());
1332 printf
1333 ("That feature is not tested further by this program.\n");
1337 Pause ();
1338 /*=============================================*/
1339 Milestone = 35;
1340 /*=============================================*/
1341 if (Radix >= Two)
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.");
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
1372 GMult = No;
1373 TstCond (Serious, false, "* lacks a Guard Digit, so 1*X != X");
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
1404 GDiv = No;
1405 TstCond (Serious, false, "Division lacks a Guard Digit, so X/1 != X");
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
1430 GAddSub = No;
1431 TstCond (Serious, false,
1432 "- lacks Guard Digit, so cancellation is obscured");
1434 if (F9 != One && F9 - One >= Zero)
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");
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;
1457 AInvrse = Radix;
1460 X = AInvrse;
1461 AInvrse = AInvrse / A1;
1463 while (!(FLOOR (AInvrse) != AInvrse));
1464 Done = (X == One) || (A1 > Three);
1465 if (!Done)
1466 A1 = Nine + One;
1468 while (!(Done));
1469 if (X == One)
1470 A1 = Radix;
1471 AInvrse = One / A1;
1472 X = A1;
1473 Y = AInvrse;
1474 Done = false;
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;
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))
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))
1515 RMult = Rounded;
1516 printf ("Multiplication appears to round correctly.\n");
1518 else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
1519 && (T < Zero) && (StickyBit + U2 == Zero) && (Y1 < Half))
1521 RMult = Chopped;
1522 printf ("Multiplication appears to chop.\n");
1524 else
1525 printf ("* is neither chopped nor correctly rounded.\n");
1526 if ((RMult == Rounded) && (GMult == No))
1527 notify ("Multiplication");
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)))
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))
1564 RDiv = Rounded;
1565 printf ("Division appears to round correctly.\n");
1566 if (GDiv == No)
1567 notify ("Division");
1569 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
1570 && (Y2 < Zero) && (Y1 - Half < F9 - Half))
1572 RDiv = Chopped;
1573 printf ("Division appears to chop.\n");
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))
1594 RAddSub = Chopped;
1595 printf ("Add/Subtract appears to be chopped.\n");
1597 if (GAddSub == Yes)
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))
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))
1615 RAddSub = Rounded;
1616 printf ("Addition/Subtraction appears to round correctly.\n");
1617 if (GAddSub == No)
1618 notify ("Add/Subtract");
1620 else
1621 printf ("Addition/Subtraction neither rounds nor chops.\n");
1623 else
1624 printf ("Addition/Subtraction neither rounds nor chops.\n");
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)
1636 S = Zero;
1637 BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n");
1639 StickyBit = Zero;
1640 if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
1641 && (RMult == Rounded) && (RDiv == Rounded)
1642 && (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2))
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))
1651 Z = T + Y;
1652 Y = Z - X;
1653 if ((Z - T >= U2) && (Y - T == Zero))
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))
1661 Z = (Half - U1) * U1;
1662 T = F9 - Z;
1663 Q = F9 - Y;
1664 if ((T - F9 == Zero) && (F9 - U1 - Q == Zero))
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)
1673 if (Radix != Two)
1675 X = Two + U2;
1676 Y = X / Two;
1677 if ((Y - One == Zero))
1678 StickyBit = S;
1680 else
1681 StickyBit = S;
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;
1707 X = Random ();
1708 Y = Random ();
1709 Z9 = Y * X;
1710 Z = X * Y;
1711 Z9 = Z - Z9;
1712 I = I + 1;
1714 while (!((I > NoTrials) || (Z9 != Zero)));
1715 if (I == NoTrials)
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);
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))
1757 X = Y;
1758 Y = Radix * Y;
1760 while (!((Y - X >= NoTrials)));
1761 OneUlp = X * U2;
1762 I = 1;
1763 while (I <= NoTrials)
1765 X = X + One;
1766 SqXMinX (Defect);
1767 if (J > Zero)
1768 break;
1769 I = I + 1;
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))
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
1788 Q = FLOOR (Q + Half);
1789 if (!(I > 0 || Radix == Q * Q))
1790 Monot = true;
1791 else if (I > 0)
1793 if (I > 1)
1794 Monot = true;
1795 else
1797 Y = Y * BInvrse;
1798 X = Y - U1;
1799 Z = Y + U1;
1802 else
1804 Y = Q;
1805 X = Y - U2;
1806 Z = Y + U2;
1810 if (Monot)
1811 printf ("sqrt has passed a test for Monotonicity.\n");
1812 else
1814 BadCond (Defect, "");
1815 printf ("sqrt(X) is non-monotonic for X near %s .\n", Y.str());
1817 /*=============================================*/
1818 Milestone = 110;
1819 /*=============================================*/
1820 printf ("Seeking Underflow thresholds UfThold and E0.\n");
1821 D = U1;
1822 if (Precision != FLOOR (Precision))
1824 D = BInvrse;
1825 X = Precision;
1828 D = D * BInvrse;
1829 X = X - One;
1831 while (X > Zero);
1833 Y = One;
1834 Z = D;
1835 /* ... D is power of 1/Radix < 1. */
1838 C = Y;
1839 Y = Z;
1840 Z = Y * Y;
1842 while ((Y > Z) && (Z + Z > Z));
1843 Y = C;
1844 Z = Y * D;
1847 C = Y;
1848 Y = Z;
1849 Z = Y * D;
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 */
1864 Y = E0;
1865 E0 = Z;
1866 Z = E0 * H;
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)
1877 E9 = Radix * U2;
1878 S = One + E9;
1879 D = C * S;
1880 if (D <= C)
1882 BadCond (Failure,
1883 "multiplication gets too many last digits wrong.\n");
1884 Underflow = E0;
1885 Y1 = Zero;
1886 PseudoZero = Z;
1887 Pause ();
1890 else
1892 Underflow = D;
1893 PseudoZero = Underflow * H;
1894 UfThold = Zero;
1897 Y1 = Underflow;
1898 Underflow = PseudoZero;
1899 if (E1 + E1 <= E1)
1901 Y2 = Underflow * HInvrse;
1902 E1 = FABS (Y1 - Y2);
1903 Q = Y1;
1904 if ((UfThold == Zero) && (Y1 != Y2))
1905 UfThold = Y1;
1907 PseudoZero = PseudoZero * H;
1909 while ((Underflow > PseudoZero)
1910 && (PseudoZero + PseudoZero > PseudoZero));
1912 /* Comment line 4530 .. 4560 */
1913 if (PseudoZero != Zero)
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)
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)
1928 printf ("But -PseudoZero, which should be\n");
1929 printf ("positive, isn't; it prints out as %s .\n", X.str());
1932 else
1934 BadCond (Flaw, "Underflow can stick at an allegedly positive\n");
1935 printf ("value PseudoZero that prints out as %s .\n",
1936 PseudoZero.str());
1938 TstPtUf ();
1940 /*=============================================*/
1941 Milestone = 120;
1942 /*=============================================*/
1943 if (CInvrse * Y > CInvrse * Y1)
1945 S = H * S;
1946 E0 = Underflow;
1948 if (!((E1 == Zero) || (E1 == E0)))
1950 BadCond (Defect, "");
1951 if (E1 < E0)
1953 printf ("Products underflow at a higher");
1954 printf (" threshold than differences.\n");
1955 if (PseudoZero == Zero)
1956 E0 = E1;
1958 else
1960 printf ("Difference underflows at a higher");
1961 printf (" threshold than products.\n");
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)
1978 case 1:
1979 UfThold = Underflow;
1980 if ((CInvrse * Q) != ((CInvrse * Y) * S))
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");
1989 Pause ();
1990 break;
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;
2001 case 3:
2002 X = X;
2003 break;
2005 case 4:
2006 if ((Q == UfThold) && (E1 == E0) && (FABS (UfThold - E1 / E9) <= E1))
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);
2018 if (UfNGrad)
2020 printf ("\n");
2021 if (setjmp (ovfl_buf))
2023 printf ("Underflow / UfThold failed!\n");
2024 R = H + H;
2026 else
2027 R = SQRT (Underflow / UfThold);
2028 if (R <= H)
2030 Z = R * UfThold;
2031 X = Z * (One + R * H * (One + H));
2033 else
2035 Z = UfThold;
2036 X = Z * (One + H * H * (One + H));
2038 if (!((X == Z) || (X - Z != Zero)))
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());
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)
2064 if (Y > E0)
2066 BadCond (Defect, "");
2067 I = 5;
2069 else
2071 BadCond (Serious, "");
2072 I = 4;
2074 printf ("Range is too narrow; U1^%d Underflows.\n", I);
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))
2088 BadCond (Serious, "trap on underflow.\n");
2090 else
2092 V9 = POW (HInvrse, Y2);
2093 printf (" %s .\n", V9.str());
2094 if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold)))
2096 BadCond (Serious, "this is not between 0 and underflow\n");
2097 printf (" threshold = %s .\n", UfThold.str());
2099 else if (!(V9 > UfThold * (One + E9)))
2100 printf ("This computed value is O.K.\n");
2101 else
2103 BadCond (Defect, "this is not between 0 and underflow\n");
2104 printf (" threshold = %s .\n", UfThold.str());
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))
2117 I = 0;
2118 V9 = Y;
2119 goto overflow;
2123 V = Y;
2124 Y = V9;
2125 V9 = HInvrse * Y;
2127 while (V9 < Y);
2128 I = 1;
2129 overflow:
2130 Z = V9;
2131 printf ("Can `Z = -Y' overflow?\n");
2132 printf ("Trying it on Y = %s .\n", Y.str());
2133 V9 = -Y;
2134 V0 = V9;
2135 if (V - Y == V + V0)
2136 printf ("Seems O.K.\n");
2137 else
2139 printf ("finds a ");
2140 BadCond (Flaw, "-(-Y) differs from Y.\n");
2142 if (Z != Y)
2144 BadCond (Serious, "");
2145 printf ("overflow past %s\n\tshrinks to %s .\n", Y.str(), Z.str());
2147 if (I)
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;
2158 else
2160 V = Y * (HInvrse * U2 - HInvrse);
2161 V = V + ((One - HInvrse) * U2) * Y;
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))
2180 BadCond (Failure, "Comparisons involving ");
2181 printf ("+-%s, +-%s\nand +-%s are confused by Overflow.",
2182 V.str(), V0.str(), UfThold.str());
2184 /*=============================================*/
2185 Milestone = 175;
2186 /*=============================================*/
2187 printf ("\n");
2188 for (Indx = 1; Indx <= 3; ++Indx)
2190 switch (Indx)
2192 case 1:
2193 Z = UfThold;
2194 break;
2195 case 2:
2196 Z = E0;
2197 break;
2198 case 3:
2199 Z = PseudoZero;
2200 break;
2202 if (Z != Zero)
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());
2218 /*=============================================*/
2219 Milestone = 180;
2220 /*=============================================*/
2221 for (Indx = 1; Indx <= 2; ++Indx)
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)))
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());
2241 /*=============================================*/
2242 Milestone = 190;
2243 /*=============================================*/
2244 Pause ();
2245 X = UfThold * V;
2246 Y = Radix * Radix;
2247 if (X * Y < One || X > Y)
2249 if (X * Y < U1 || X > Y / U1)
2250 BadCond (Defect, "Badly");
2251 else
2252 BadCond (Flaw, "");
2254 printf (" unbalanced range; UfThold * V = %s\n\t%s\n",
2255 X.str(), "is too far from 1.\n");
2257 /*=============================================*/
2258 Milestone = 200;
2259 /*=============================================*/
2260 for (Indx = 1; Indx <= 5; ++Indx)
2262 X = F9;
2263 switch (Indx)
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;
2277 Y = X;
2278 if (setjmp (ovfl_buf))
2279 printf (" X / X traps when X = %s\n", X.str());
2280 else
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());
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");
2311 static const char *msg[] = {
2312 "FAILUREs encountered =",
2313 "SERIOUS DEFECTs discovered =",
2314 "DEFECTs discovered =",
2315 "FLAWs discovered ="
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]);
2322 printf ("\n");
2323 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] + ErrCnt[Flaw]) > 0)
2325 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] == 0)
2326 && (ErrCnt[Flaw] > 0))
2328 printf ("The arithmetic diagnosed seems ");
2329 printf ("Satisfactory though flawed.\n");
2331 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0) && (ErrCnt[Defect] > 0))
2333 printf ("The arithmetic diagnosed may be Acceptable\n");
2334 printf ("despite inconvenient Defects.\n");
2336 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0)
2338 printf ("The arithmetic diagnosed has ");
2339 printf ("unacceptable Serious Defects.\n");
2341 if (ErrCnt[Failure] > 0)
2343 printf ("Potentially fatal FAILURE may have spoiled this");
2344 printf (" program's subsequent diagnoses.\n");
2347 else
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
2355 if (StickyBit >= One &&
2356 (Radix - Two) * (Radix - Nine - One) == Zero)
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
2370 printf (",\nexcept for possibly Double Rounding");
2371 printf (" during Gradual Underflow.\n");
2374 printf ("The arithmetic diagnosed appears to be Excellent!\n");
2377 printf ("END OF TEST.\n");
2378 return 0;
2381 template<typename FLOAT>
2382 FLOAT
2383 Paranoia<FLOAT>::Sign (FLOAT X)
2385 return X >= FLOAT (long (0)) ? 1 : -1;
2388 template<typename FLOAT>
2389 void
2390 Paranoia<FLOAT>::Pause ()
2392 if (do_pause)
2394 fputs ("Press return...", stdout);
2395 fflush (stdout);
2396 getchar();
2398 printf ("\nDiagnosis resumes after milestone Number %d", Milestone);
2399 printf (" Page: %d\n\n", PageNo);
2400 ++Milestone;
2401 ++PageNo;
2404 template<typename FLOAT>
2405 void
2406 Paranoia<FLOAT>::TstCond (int K, int Valid, const char *T)
2408 if (!Valid)
2410 BadCond (K, T);
2411 printf (".\n");
2415 template<typename FLOAT>
2416 void
2417 Paranoia<FLOAT>::BadCond (int K, const char *T)
2419 static const char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
2421 ErrCnt[K] = ErrCnt[K] + 1;
2422 printf ("%s: %s", msg[K], T);
2425 /* Random computes
2426 X = (Random1 + Random9)^5
2427 Random1 = X - FLOOR(X) + 0.000005 * X;
2428 and returns the new value of Random1. */
2430 template<typename FLOAT>
2431 FLOAT
2432 Paranoia<FLOAT>::Random ()
2434 FLOAT X, Y;
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);
2445 template<typename FLOAT>
2446 void
2447 Paranoia<FLOAT>::SqXMinX (int ErrKind)
2449 FLOAT XA, XB;
2451 XB = X * BInvrse;
2452 XA = X - XB;
2453 SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp;
2454 if (SqEr != Zero)
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");
2468 template<typename FLOAT>
2469 void
2470 Paranoia<FLOAT>::NewD ()
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)
2478 Z = -Z;
2479 Z1 = -Z1;
2481 D = Radix * D;
2484 template<typename FLOAT>
2485 void
2486 Paranoia<FLOAT>::SR3750 ()
2488 if (!((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2)))
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;
2504 template<typename FLOAT>
2505 void
2506 Paranoia<FLOAT>::IsYeqX ()
2508 if (Y != X)
2510 if (N <= 0)
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());
2521 N = N + 1; /* ... count discrepancies. */
2525 template<typename FLOAT>
2526 void
2527 Paranoia<FLOAT>::PrintIfNPositive ()
2529 if (N > 0)
2530 printf ("Similar discrepancies have occurred %d times.\n", N);
2533 template<typename FLOAT>
2534 void
2535 Paranoia<FLOAT>::TstPtUf ()
2537 N = 0;
2538 if (Z != Zero)
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)
2548 printf ("This is O.K., provided Over/Underflow");
2549 printf (" has NOT just been signaled.\n");
2551 else
2553 if ((Q9 < One) || (Q9 > Two))
2555 very_serious:
2556 N = 1;
2557 ErrCnt[Serious] = ErrCnt[Serious] + 1;
2558 printf ("This is a VERY SERIOUS DEFECT!\n");
2560 else
2562 N = 1;
2563 ErrCnt[Defect] = ErrCnt[Defect] + 1;
2564 printf ("This is a DEFECT!\n");
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))
2574 if (N > 0)
2575 Pause ();
2577 else
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)
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());
2595 Pause ();
2600 template<typename FLOAT>
2601 void
2602 Paranoia<FLOAT>::notify (const char *s)
2604 printf ("%s test appears to be inconsistent...\n", s);
2605 printf (" PLEASE NOTIFY KARPINKSI!\n");
2608 /* ====================================================================== */
2610 int main(int ac, char **av)
2612 setbuf(stdout, NULL);
2613 setbuf(stderr, NULL);
2615 while (1)
2616 switch (getopt (ac, av, "pvg:fdl"))
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':
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
2649 int i, n = sizeof (fmts)/sizeof(*fmts);
2651 for (i = 0; i < n; ++i)
2652 if (strcmp (fmts[i].name, optarg) == 0)
2653 break;
2655 if (i == n)
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;
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;
2669 Paranoia<real_c_float>().main();
2670 break;
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;
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;
2695 /* GCC stuff referenced by real.o. */
2697 extern "C" void
2698 fancy_abort ()
2700 abort ();
2703 int target_flags = 0;
2705 extern "C" int
2706 floor_log2_wide (unsigned HOST_WIDE_INT x)
2708 int log = -1;
2709 while (x != 0)
2710 log++,
2711 x >>= 1;
2712 return log;