OSDN Git Service

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