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,
17 Electrical Engineering & Computer Science Dept.
18 University of California
19 Berkeley, California 94720
22 converted to Pascal by:
24 National Physical Laboratory
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
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
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
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
85 F1 MinusOne O5 Five X9 Random1
91 G2 GDiv Q9 Z0 PseudoZero
94 H1 HInverse R2 RDiv Z9
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:
117 170- 250 Instructions
119 480- 670 Characteristics
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. */
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
152 #include "ansidecl.h"
153 #include "auto-host.h"
156 #undef EXTRA_MODES_FILE
159 typedef struct rtx_def
*rtx
;
161 typedef struct rtvec_def
*rtvec
;
163 typedef union tree_node
*tree
;
165 #define DEFTREECODE(SYM, STRING, TYPE, NARGS) SYM,
168 LAST_AND_UNUSED_TREE_CODE
172 #define ENUM_BITFIELD(X) enum X
180 /* We never produce signals from the library. Thus setjmp need do nothing. */
182 #define setjmp(x) (0)
184 static bool verbose
= false;
185 static int verbose_index
= 0;
187 /* ====================================================================== */
188 /* The implementation of the abstract floating point class based on gcc's
189 real.c. I.e. the object of this excersize. Templated so that we can
195 static const enum machine_mode MODE
= SFmode
;
198 long image
[128 / 32];
200 void from_long(long);
201 void from_str(const char *);
202 void binop(int code
, const real_c_float
&);
204 bool cmp(int code
, const real_c_float
&) const;
211 real_c_float(const char *s
)
213 real_c_float(const real_c_float
&b
)
214 { memcpy(image
, b
.image
, sizeof(image
)); }
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; }
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; }
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
; }
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
); }
244 const char * str () const;
245 const char * hex () const;
246 long integer () const;
252 real_c_float::from_long (long l
)
256 real_from_integer (&f
, MODE
, l
, l
< 0 ? -1 : 0, 0);
257 real_to_target (image
, &f
, MODE
);
261 real_c_float::from_str (const char *s
)
266 if (*p
== '-' || *p
== '+')
268 if (strcasecmp(p
, "inf") == 0)
272 real_arithmetic (&f
, NEGATE_EXPR
, &f
, NULL
);
274 else if (strcasecmp(p
, "nan") == 0)
275 real_nan (&f
, "", 1, MODE
);
277 real_from_string (&f
, s
);
279 real_to_target (image
, &f
, MODE
);
283 real_c_float::binop (int code
, const real_c_float
&b
)
285 REAL_VALUE_TYPE ai
, bi
, ri
;
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
);
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
;
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);
307 symbol_for_code
= '+';
310 symbol_for_code
= '-';
313 symbol_for_code
= '*';
316 symbol_for_code
= '/';
322 fprintf (stderr
, "%6d: %s %c %s = %s\n", verbose_index
++,
323 ab
, symbol_for_code
, bb
, rb
);
328 real_c_float::unop (int code
)
330 REAL_VALUE_TYPE ai
, ri
;
332 real_from_target (&ai
, image
, MODE
);
333 real_arithmetic (&ri
, code
, &ai
, NULL
);
334 real_to_target (image
, &ri
, MODE
);
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
;
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);
350 symbol_for_code
= "-";
353 symbol_for_code
= "abs ";
359 fprintf (stderr
, "%6d: %s%s = %s\n", verbose_index
++,
360 symbol_for_code
, ab
, rb
);
365 real_c_float::cmp (int code
, const real_c_float
&b
) const
367 REAL_VALUE_TYPE ai
, bi
;
370 real_from_target (&ai
, image
, MODE
);
371 real_from_target (&bi
, b
.image
, MODE
);
372 ret
= real_compare (code
, &ai
, &bi
);
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
;
381 real_to_hexadecimal (ab
, &ai
, sizeof(ab
), digits
, 0);
382 real_to_hexadecimal (bb
, &bi
, sizeof(bb
), digits
, 0);
387 symbol_for_code
= "<";
390 symbol_for_code
= "<=";
393 symbol_for_code
= "==";
396 symbol_for_code
= "!=";
399 symbol_for_code
= ">=";
402 symbol_for_code
= ">";
408 fprintf (stderr
, "%6d: %s %s %s = %s\n", verbose_index
++,
409 ab
, symbol_for_code
, bb
, (ret
? "true" : "false"));
416 real_c_float::str() const
419 const real_format
*fmt
= real_format_for_mode
[MODE
- QFmode
];
420 const int digits
= int(fmt
->p
* fmt
->log2_b
* .30102999566398119521 + 1);
422 real_from_target (&f
, image
, MODE
);
423 char *buf
= new char[digits
+ 10];
424 real_to_decimal (buf
, &f
, digits
+10, digits
, 0);
430 real_c_float::hex() const
433 const real_format
*fmt
= real_format_for_mode
[MODE
- QFmode
];
434 const int digits
= (fmt
->p
* fmt
->log2_b
+ 3) / 4;
436 real_from_target (&f
, image
, MODE
);
437 char *buf
= new char[digits
+ 10];
438 real_to_hexadecimal (buf
, &f
, digits
+10, digits
, 0);
444 real_c_float::integer() const
447 real_from_target (&f
, image
, MODE
);
448 return real_to_integer (&f
);
452 real_c_float::exp() const
455 real_from_target (&f
, image
, MODE
);
456 return real_exponent (&f
);
460 real_c_float::ldexp (int exp
)
464 real_from_target (&ai
, image
, MODE
);
465 real_ldexp (&ai
, &ai
, exp
);
466 real_to_target (image
, &ai
, MODE
);
469 /* ====================================================================== */
470 /* An implementation of the abstract floating point class that uses native
471 arithmetic. Exists for reference and debugging. */
477 // Force intermediate results back to memory.
480 static T
from_str (const char *);
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);
491 native_float(const char *s
)
492 { image
= from_str(s
); }
493 native_float(const native_float
&b
)
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; }
503 const native_float
& operator+= (const native_float
&b
)
505 image
= verbose_binop(image
, '+', b
.image
, image
+ b
.image
);
508 const native_float
& operator-= (const native_float
&b
)
510 image
= verbose_binop(image
, '-', b
.image
, image
- b
.image
);
513 const native_float
& operator*= (const native_float
&b
)
515 image
= verbose_binop(image
, '*', b
.image
, image
* b
.image
);
518 const native_float
& operator/= (const native_float
&b
)
520 image
= verbose_binop(image
, '/', b
.image
, image
/ b
.image
);
524 native_float
operator- () const
527 r
.image
= verbose_unop("-", image
, -image
);
530 native_float
abs () const
533 r
.image
= verbose_unop("abs ", image
, do_abs(image
));
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
); }
550 const char * str () const;
551 const char * hex () const;
552 long integer () const
553 { return long(image
); }
560 native_float
<T
>::from_str (const char *s
)
562 return strtold (s
, NULL
);
567 native_float
<float>::from_str (const char *s
)
569 return strtof (s
, NULL
);
574 native_float
<double>::from_str (const char *s
)
576 return strtod (s
, NULL
);
581 native_float
<T
>::do_abs (T image
)
583 return fabsl (image
);
588 native_float
<float>::do_abs (float image
)
590 return fabsf (image
);
595 native_float
<double>::do_abs (double image
)
602 native_float
<T
>::verbose_binop (T a
, char symbol
, T b
, T r
)
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
);
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
);
622 native_float
<T
>::verbose_unop (const char *symbol
, T a
, T r
)
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
);
631 fprintf (stderr
, "%6d: %s%.*La = %.*La\n", verbose_index
++,
632 symbol
, digits
, (long double)a
, digits
, (long double)r
);
640 native_float
<T
>::verbose_cmp (T a
, const char *symbol
, T b
, bool r
)
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"));
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"));
660 native_float
<T
>::str() const
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
);
667 sprintf (buf
, "%.*Le", digits
- 1, (long double) image
);
674 native_float
<T
>::hex() const
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
);
681 sprintf (buf
, "%.*La", digits
- 1, (long double) image
);
688 native_float
<T
>::exp() const
697 native_float
<T
>::ldexp (int exp
)
699 image
= ldexpl (image
, exp
);
704 native_float
<float>::ldexp (int exp
)
706 image
= ldexpf (image
, exp
);
711 native_float
<double>::ldexp (int exp
)
713 image
= ::ldexp (image
, exp
);
716 /* ====================================================================== */
717 /* Some libm routines that Paranoia expects to be available. */
719 template<typename FLOAT
>
721 FABS (const FLOAT
&f
)
726 template<typename FLOAT
, typename RHS
>
728 operator+ (const FLOAT
&a
, const RHS
&b
)
730 return FLOAT(a
) += FLOAT(b
);
733 template<typename FLOAT
, typename RHS
>
735 operator- (const FLOAT
&a
, const RHS
&b
)
737 return FLOAT(a
) -= FLOAT(b
);
740 template<typename FLOAT
, typename RHS
>
742 operator* (const FLOAT
&a
, const RHS
&b
)
744 return FLOAT(a
) *= FLOAT(b
);
747 template<typename FLOAT
, typename RHS
>
749 operator/ (const FLOAT
&a
, const RHS
&b
)
751 return FLOAT(a
) /= FLOAT(b
);
754 template<typename FLOAT
>
756 FLOOR (const FLOAT
&f
)
758 /* ??? This is only correct when F is representable as an integer. */
759 long i
= f
.integer();
769 template<typename FLOAT
>
771 SQRT (const FLOAT
&f
)
774 FLOAT zero
= long(0);
788 z
.ldexp (-f
.exp() / 2);
790 diff2
= FABS (z
* z
- f
);
794 t
= (f
/ (two
* z
)) + (z
/ two
);
795 diff
= FABS (t
* t
- f
);
803 #elif defined(NO_LONG_DOUBLE)
807 d
= strtod (f
.hex(), NULL
);
809 sprintf(buf
, "%.35a", d
);
816 ld
= strtold (f
.hex(), NULL
);
818 sprintf(buf
, "%.35La", ld
);
824 template<typename FLOAT
>
829 FLOAT zero
= long(0);
837 int exp
= x
.exp() - 1;
848 FLOAT term
= y
/ FLOAT (n
);
849 FLOAT next
= sum
+ term
;
858 sum
+= FLOAT (exp
) * FLOAT(".69314718055994530941");
861 #elif defined (NO_LONG_DOUBLE)
865 d
= strtod (x
.hex(), NULL
);
867 sprintf(buf
, "%.35a", d
);
874 ld
= strtold (x
.hex(), NULL
);
876 sprintf(buf
, "%.35La", ld
);
882 template<typename FLOAT
>
887 #ifdef NO_LONG_DOUBLE
891 d
= strtod (x
.hex(), NULL
);
893 sprintf(buf
, "%.35a", d
);
900 ld
= strtold (x
.hex(), NULL
);
902 sprintf(buf
, "%.35La", ld
);
908 template<typename FLOAT
>
910 POW (const FLOAT
&base
, const FLOAT
&exp
)
913 #ifdef NO_LONG_DOUBLE
917 d1
= strtod (base
.hex(), NULL
);
918 d2
= strtod (exp
.hex(), NULL
);
920 sprintf(buf
, "%.35a", d1
);
924 long double ld1
, ld2
;
927 ld1
= strtold (base
.hex(), NULL
);
928 ld2
= strtold (exp
.hex(), NULL
);
929 ld1
= powl (ld1
, ld2
);
930 sprintf(buf
, "%.35La", ld1
);
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. */
940 int NoTrials
= 20; /*Number of tests for commutativity. */
941 bool do_pause
= false;
943 enum Guard
{ No
, Yes
};
944 enum Rounding
{ Other
, Rounded
, Chopped
};
945 enum Class
{ Failure
, Serious
, Defect
, Flaw
};
947 template<typename FLOAT
>
950 FLOAT Radix
, BInvrse
, RadixD2
, BMinusU2
;
952 /* Small floating point constants. */
968 /* Declarations of Variables. */
974 FLOAT E0
, E1
, Exp2
, E3
, MinSqEr
;
975 FLOAT SqEr
, MaxSqEr
, E9
;
985 FLOAT T
, Underflow
, S
;
986 FLOAT OneUlp
, UfThold
, U1
, U2
;
989 FLOAT X
, X1
, X2
, X8
, Random1
;
990 FLOAT Y
, Y1
, Y2
, Random2
;
991 FLOAT Z
, PseudoZero
, Z1
, Z2
, Z9
;
996 Guard GMult
, GDiv
, GAddSub
;
997 Rounding RMult
, RDiv
, RAddSub
, RSqrt
;
998 int Break
, Done
, NotMonot
, Monot
, Anomaly
, IEEE
, SqRWrng
, UfNGrad
;
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 */
1009 void BadCond (int, const char *);
1011 void TstCond (int, int, const char *);
1012 void notify (const char *);
1015 void PrintIfNPositive ();
1019 // Pretend we're bss.
1020 Paranoia() { memset(this, 0, sizeof (*this)); }
1023 template<typename FLOAT
>
1025 Paranoia
<FLOAT
>::main()
1027 /* First two assignments use integer right-hand sides. */
1036 TwentySeven
= long(27);
1037 ThirtyTwo
= long(32);
1038 TwoForty
= long(240);
1039 MinusOne
= long(-1);
1041 OneAndHalf
= "0x3p-1";
1042 ErrCnt
[Failure
] = 0;
1043 ErrCnt
[Serious
] = 0;
1047 /*=============================================*/
1049 /*=============================================*/
1050 printf ("Program is now RUNNING tests on small integers:\n");
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");
1060 ErrCnt
[Failure
] = ErrCnt
[Failure
] + 1;
1061 printf ("Comparison alleges that -0.0 is Non-zero!\n");
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");
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");
1079 TstCond (Failure
, Half
+ MinusOne
+ Half
== Zero
, "1/2 + (-1) + 1/2 != 0");
1081 /*=============================================*/
1083 /*=============================================*/
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
),
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");
1101 if (ErrCnt
[Failure
] == 0)
1103 printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
1106 printf ("Searching for Radix and Precision.\n");
1115 while (MinusOne
+ FABS (Y
) < Zero
);
1116 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ... */
1125 while (Radix
== Zero
);
1128 printf ("Radix = %s .\n", Radix
.str());
1134 Precision
= Precision
+ One
;
1138 while ((Y
- W
) == One
);
1140 /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
1144 printf ("Closest relative separation found is U1 = %s .\n\n", U1
.str());
1145 printf ("Recalculating radix and precision\n ");
1147 /*save old values */
1157 X
= FABS (X
- Third
);
1161 /*... now X = (unknown no.) ulps of 1+... */
1165 Y
= Half
* U2
+ ThirtyTwo
* U2
* U2
;
1169 while (!((U2
<= X
) || (X
<= Zero
)));
1171 /*... now U2 == 1 ulp of 1 + ... */
1180 /*... now X == (unknown no.) ulps of 1 -... */
1184 Y
= Half
* U1
+ ThirtyTwo
* U1
* U1
;
1190 while (!((U1
<= X
) || (X
<= Zero
)));
1191 /*... now U1 == 1 ulp of 1 - ... */
1193 printf ("confirms closest relative separation U1 .\n");
1195 printf ("gets better closest relative separation U1 = %s .\n", U1
.str());
1197 F9
= (Half
- U1
) + Half
;
1199 Radix
= FLOOR (FLOAT ("0.01") + U2
/ U1
);
1201 printf ("Radix confirmed.\n");
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 /*=============================================*/
1210 /*=============================================*/
1211 TstCond (Failure
, F9
- Half
< Half
,
1212 "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
1217 TstCond (Failure
, (X
!= One
)
1218 || (Z
== Zero
), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
1221 /*=============================================*/
1223 /*=============================================*/
1224 /*... BMinusU2 = nextafter(Radix, 0) */
1225 BMinusU2
= Radix
- One
;
1226 BMinusU2
= (BMinusU2
- U2
) + One
;
1227 /* Purify Integers */
1230 X
= -TwoForty
* LOG (U1
) / LOG (Radix
);
1231 Y
= FLOOR (Half
+ X
);
1232 if (FABS (X
- Y
) * Four
< One
)
1234 Precision
= X
/ TwoForty
;
1235 Y
= FLOOR (Half
+ Precision
);
1236 if (FABS (Precision
- Y
) * TwoForty
< Half
)
1239 if ((Precision
!= FLOOR (Precision
)) || (Radix
== One
))
1241 printf ("Precision cannot be characterized by an Integer number\n");
1243 ("of significant digits but, by itself, this is a minor flaw.\n");
1247 ("logarithmic encoding has precision characterized solely by U1.\n");
1249 printf ("The number of significant digits of the Radix is %s .\n",
1251 TstCond (Serious
, U2
* Nine
* Nine
* TwoForty
< One
,
1252 "Precision worse than 5 decimal figures ");
1253 /*=============================================*/
1255 /*=============================================*/
1256 /* Test for extra-precise subexpressions */
1257 X
= FABS (((Four
/ Three
- One
) - One
/ Four
) * Three
- One
/ Four
);
1261 X
= (One
+ (Half
* Z2
+ ThirtyTwo
* Z2
* Z2
)) - One
;
1263 while (!((Z2
<= X
) || (X
<= Zero
)));
1264 X
= Y
= Z
= FABS ((Three
/ Four
- Two
/ Three
) * Three
- One
/ Four
);
1268 Z
= (One
/ Two
- ((One
/ Two
- (Half
* Z1
+ ThirtyTwo
* Z1
* Z1
))
1269 + One
/ Two
)) + One
/ Two
;
1271 while (!((Z1
<= Z
) || (Z
<= Zero
)));
1278 (Half
- ((Half
- (Half
* Y1
+ ThirtyTwo
* Y1
* Y1
)) + Half
)) +
1281 while (!((Y1
<= Y
) || (Y
<= Zero
)));
1283 X
= ((Half
* X1
+ ThirtyTwo
* X1
* X1
) - F9
) + F9
;
1285 while (!((X1
<= X
) || (X
<= Zero
)));
1286 if ((X1
!= Y1
) || (X1
!= Z1
))
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");
1298 if ((Z1
!= U1
) || (Z2
!= U2
))
1300 if ((Z1
>= U1
) || (Z2
>= U2
))
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());
1309 if ((Z1
<= Zero
) || (Z2
<= Zero
))
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");
1316 if (Z1
!= Z2
|| Z1
> Zero
)
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());
1331 ("That feature is not tested further by this program.\n");
1336 /*=============================================*/
1338 /*=============================================*/
1341 X
= W
/ (Radix
* Radix
);
1346 TstCond (Failure
, X
== U2
,
1347 "Subtraction is not normalized X=Y,X+Z != Y+Z!");
1349 printf ("Subtraction appears to be normalized, as it should be.");
1351 printf ("\nChecking for guard digit in *, /, and -.\n");
1364 X
= X
* (Radix
- One
);
1365 T
= T
* (Radix
- One
);
1366 if ((X
== Zero
) && (Y
== Zero
) && (Z
== Zero
) && (T
== Zero
))
1371 TstCond (Serious
, false, "* lacks a Guard Digit, so 1*X != X");
1375 Y
= FABS ((X
+ Z
) - X
* X
) - U2
;
1377 Z
= FABS ((X
- U2
) - X
* X
) - U1
;
1378 TstCond (Failure
, (Y
<= Zero
)
1379 && (Z
<= Zero
), "* gets too many final digits wrong.\n");
1387 T
= Nine
/ TwentySeven
;
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");
1398 if ((X
== Zero
) && (Y
== Zero
) && (Z
== Zero
))
1403 TstCond (Serious
, false, "Division lacks a Guard Digit, so X/1 != X");
1405 X
= One
/ (One
+ U2
);
1406 Y
= X
- Half
- Half
;
1407 TstCond (Serious
, Y
< Zero
, "Computed value of 1/1.000..1 >= 1");
1409 Y
= One
+ Radix
* U2
;
1413 StickyBit
= T
/ Radix
;
1416 TstCond (Failure
, X
== Zero
&& Y
== Zero
,
1417 "* and/or / gets too many last digits wrong");
1422 Z
= Radix
- BMinusU2
;
1424 if ((X
== U1
) && (Y
== U1
) && (Z
== U2
) && (T
== U2
))
1429 TstCond (Serious
, false,
1430 "- lacks Guard Digit, so cancellation is obscured");
1432 if (F9
!= One
&& F9
- One
>= Zero
)
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");
1439 if (GMult
== Yes
&& GDiv
== Yes
&& GAddSub
== Yes
)
1441 (" *, /, and - appear to have guard digits, as they should.\n");
1442 /*=============================================*/
1444 /*=============================================*/
1446 printf ("Checking rounding on multiply, divide and add/subtract.\n");
1450 RadixD2
= Radix
/ Two
;
1459 AInvrse
= AInvrse
/ A1
;
1461 while (!(FLOOR (AInvrse
) != AInvrse
));
1462 Done
= (X
== One
) || (A1
> Three
);
1476 TstCond (Failure
, Z
== Half
, "X * (1/X) differs from 1");
1484 X
= OneAndHalf
- U2
;
1485 Y
= OneAndHalf
+ U2
;
1494 if ((X
== Zero
) && (Y
== Zero
) && (Z
== Zero
) && (T
<= Zero
))
1496 X
= (OneAndHalf
+ U2
) * Y2
;
1497 Y
= OneAndHalf
- U2
- U2
;
1498 Z
= OneAndHalf
+ U2
+ U2
;
1499 T
= (OneAndHalf
- U2
) * Y1
;
1504 Y
= (U2
- Y
) + StickyBit
;
1505 Z
= S
- (Z
+ U2
+ U2
);
1506 StickyBit
= (Y2
+ U2
) * Y1
;
1508 StickyBit
= StickyBit
- Y2
;
1510 if ((X
== Zero
) && (Y
== Zero
) && (Z
== Zero
) && (T
== Zero
)
1511 && (StickyBit
== Zero
) && (Y1
== Half
))
1514 printf ("Multiplication appears to round correctly.\n");
1516 else if ((X
+ U2
== Zero
) && (Y
< Zero
) && (Z
+ U2
== Zero
)
1517 && (T
< Zero
) && (StickyBit
+ U2
== Zero
) && (Y1
< Half
))
1520 printf ("Multiplication appears to chop.\n");
1523 printf ("* is neither chopped nor correctly rounded.\n");
1524 if ((RMult
== Rounded
) && (GMult
== No
))
1525 notify ("Multiplication");
1528 printf ("* is neither chopped nor correctly rounded.\n");
1529 /*=============================================*/
1531 /*=============================================*/
1534 Z
= OneAndHalf
+ U2
+ U2
;
1536 T
= OneAndHalf
- U2
- U2
;
1542 Z
= Z
- (OneAndHalf
+ U2
);
1543 T
= (U2
- OneAndHalf
) + T
;
1544 if (!((X
> Zero
) || (Y
> Zero
) || (Z
> Zero
) || (T
> Zero
)))
1546 X
= OneAndHalf
/ Y2
;
1547 Y
= OneAndHalf
- U2
;
1548 Z
= OneAndHalf
+ U2
;
1550 T
= OneAndHalf
/ Y1
;
1555 Y1
= (Y2
+ U2
) / 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
))
1563 printf ("Division appears to round correctly.\n");
1565 notify ("Division");
1567 else if ((X
< Zero
) && (Y
< Zero
) && (Z
< Zero
) && (T
< Zero
)
1568 && (Y2
< Zero
) && (Y1
- Half
< F9
- Half
))
1571 printf ("Division appears to chop.\n");
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 /*=============================================*/
1581 /*=============================================*/
1582 TstCond (Failure
, ((F9
+ U1
) - Half
== Half
)
1583 && ((BMinusU2
+ U2
) - One
== Radix
- One
),
1584 "Incomplete carry-propagation in Addition");
1586 Y
= One
+ U2
* (One
- U2
);
1590 if ((X
== Zero
) && (Y
== Zero
))
1593 printf ("Add/Subtract appears to be chopped.\n");
1597 X
= (Half
+ U2
) * U2
;
1598 Y
= (Half
- U2
) * U2
;
1603 if ((X
== Zero
) && (Y
== Zero
))
1605 X
= (Half
+ U2
) * U1
;
1606 Y
= (Half
- U2
) * U1
;
1611 if ((X
== Zero
) && (Y
== Zero
))
1614 printf ("Addition/Subtraction appears to round correctly.\n");
1616 notify ("Add/Subtract");
1619 printf ("Addition/Subtraction neither rounds nor chops.\n");
1622 printf ("Addition/Subtraction neither rounds nor chops.\n");
1625 printf ("Addition/Subtraction neither rounds nor chops.\n");
1627 X
= One
+ Half
* (One
+ Half
);
1628 Y
= (One
+ U2
) * Half
;
1632 if (StickyBit
!= Zero
)
1635 BadCond (Flaw
, "(X - Y) + (Y - X) is non zero!\n");
1638 if ((GMult
== Yes
) && (GDiv
== Yes
) && (GAddSub
== Yes
)
1639 && (RMult
== Rounded
) && (RDiv
== Rounded
)
1640 && (RAddSub
== Rounded
) && (FLOOR (RadixD2
) == RadixD2
))
1642 printf ("Checking for sticky bit.\n");
1643 X
= (Half
+ U1
) * U2
;
1647 if ((Z
- One
<= Zero
) && (T
- One
>= U2
))
1651 if ((Z
- T
>= U2
) && (Y
- T
== Zero
))
1653 X
= (Half
+ U1
) * U1
;
1657 if ((Z
- One
== Zero
) && (T
- F9
== Zero
))
1659 Z
= (Half
- U1
) * U1
;
1662 if ((T
- F9
== Zero
) && (F9
- U1
- Q
== Zero
))
1664 Z
= (One
+ U2
) * OneAndHalf
;
1665 T
= (OneAndHalf
+ U2
) - Z
+ U2
;
1666 X
= One
+ Half
/ Radix
;
1667 Y
= One
+ Radix
* U2
;
1669 if (T
== Zero
&& X
+ Radix
* U2
- Z
== Zero
)
1675 if ((Y
- One
== Zero
))
1686 if (StickyBit
== One
)
1687 printf ("Sticky bit apparently used correctly.\n");
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 /*=============================================*/
1696 /*=============================================*/
1698 printf ("Does Multiplication commute? ");
1699 printf ("Testing on %d random pairs.\n", NoTrials
);
1700 Random9
= SQRT (FLOAT (3));
1712 while (!((I
> NoTrials
) || (Z9
!= Zero
)));
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
) +
1723 if (!((I
== NoTrials
) || (Z9
== Zero
)))
1724 BadCond (Defect
, "X * Y == Y * X trial fails.\n");
1726 printf (" No failures found in %d integer pairs.\n", NoTrials
);
1727 /*=============================================*/
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");
1741 OneUlp
= BInvrse
* U1
;
1748 printf ("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials
);
1758 while (!((Y
- X
>= NoTrials
)));
1761 while (I
<= NoTrials
)
1769 printf ("Test for sqrt monotonicity.\n");
1773 Z
= Radix
+ Radix
* U2
;
1776 while (!(NotMonot
|| Monot
))
1782 if ((X
> Q
) || (Q
> Z
))
1786 Q
= FLOOR (Q
+ Half
);
1787 if (!(I
> 0 || Radix
== Q
* Q
))
1809 printf ("sqrt has passed a test for Monotonicity.\n");
1812 BadCond (Defect
, "");
1813 printf ("sqrt(X) is non-monotonic for X near %s .\n", Y
.str());
1815 /*=============================================*/
1817 /*=============================================*/
1818 printf ("Seeking Underflow thresholds UfThold and E0.\n");
1820 if (Precision
!= FLOOR (Precision
))
1833 /* ... D is power of 1/Radix < 1. */
1840 while ((Y
> Z
) && (Z
+ Z
> Z
));
1849 while ((Y
> Z
) && (Z
+ Z
> Z
));
1855 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1859 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1866 while ((E0
> Z
) && (Z
+ Z
> Z
));
1881 "multiplication gets too many last digits wrong.\n");
1891 PseudoZero
= Underflow
* H
;
1896 Underflow
= PseudoZero
;
1899 Y2
= Underflow
* HInvrse
;
1900 E1
= FABS (Y1
- Y2
);
1902 if ((UfThold
== Zero
) && (Y1
!= Y2
))
1905 PseudoZero
= PseudoZero
* H
;
1907 while ((Underflow
> PseudoZero
)
1908 && (PseudoZero
+ PseudoZero
> PseudoZero
));
1910 /* Comment line 4530 .. 4560 */
1911 if (PseudoZero
!= Zero
)
1915 /* ... Test PseudoZero for "phoney- zero" violates */
1916 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1918 if (PseudoZero
<= Zero
)
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());
1926 printf ("But -PseudoZero, which should be\n");
1927 printf ("positive, isn't; it prints out as %s .\n", X
.str());
1932 BadCond (Flaw
, "Underflow can stick at an allegedly positive\n");
1933 printf ("value PseudoZero that prints out as %s .\n",
1938 /*=============================================*/
1940 /*=============================================*/
1941 if (CInvrse
* Y
> CInvrse
* Y1
)
1946 if (!((E1
== Zero
) || (E1
== E0
)))
1948 BadCond (Defect
, "");
1951 printf ("Products underflow at a higher");
1952 printf (" threshold than differences.\n");
1953 if (PseudoZero
== Zero
)
1958 printf ("Difference underflows at a higher");
1959 printf (" threshold than products.\n");
1962 printf ("Smallest strictly positive number found is E0 = %s .\n", E0
.str());
1971 if (UfThold
== Zero
)
1977 UfThold
= Underflow
;
1978 if ((CInvrse
* Q
) != ((CInvrse
* Y
) * S
))
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());
1985 (" or else multiplication gets too many last digits wrong.\n");
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());
2004 if ((Q
== UfThold
) && (E1
== E0
) && (FABS (UfThold
- E1
/ E9
) <= E1
))
2007 printf ("Underflow is gradual; it incurs Absolute Error =\n");
2008 printf ("(roundoff in UfThold) < E0.\n");
2010 Y
= Y
* (OneAndHalf
+ U2
);
2011 X
= CInvrse
* (One
+ U2
);
2019 if (setjmp (ovfl_buf
))
2021 printf ("Underflow / UfThold failed!\n");
2025 R
= SQRT (Underflow
/ UfThold
);
2029 X
= Z
* (One
+ R
* H
* (One
+ H
));
2034 X
= Z
* (One
+ H
* H
* (One
+ H
));
2036 if (!((X
== Z
) || (X
- Z
!= Zero
)))
2039 printf ("X = %s\n\tis not equal to Z = %s .\n", X
.str(), Z
.str());
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");
2051 printf ("X / Z = 1 + %s .\n", ((X
/ Z
- Half
) - Half
).str());
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");
2064 BadCond (Defect
, "");
2069 BadCond (Serious
, "");
2072 printf ("Range is too narrow; U1^%d Underflows.\n", I
);
2074 /*=============================================*/
2076 /*=============================================*/
2077 Y
= -FLOOR (Half
- TwoForty
* LOG (UfThold
) / LOG (HInvrse
)) / TwoForty
;
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
))
2086 BadCond (Serious
, "trap on underflow.\n");
2090 V9
= POW (HInvrse
, Y2
);
2091 printf (" %s .\n", V9
.str());
2092 if (!((V9
>= Zero
) && (V9
<= (Radix
+ Radix
+ E9
) * UfThold
)))
2094 BadCond (Serious
, "this is not between 0 and underflow\n");
2095 printf (" threshold = %s .\n", UfThold
.str());
2097 else if (!(V9
> UfThold
* (One
+ E9
)))
2098 printf ("This computed value is O.K.\n");
2101 BadCond (Defect
, "this is not between 0 and underflow\n");
2102 printf (" threshold = %s .\n", UfThold
.str());
2105 /*=============================================*/
2107 /*=============================================*/
2109 printf ("Searching for Overflow threshold:\n");
2110 printf ("This may generate an error.\n");
2113 if (setjmp (ovfl_buf
))
2129 printf ("Can `Z = -Y' overflow?\n");
2130 printf ("Trying it on Y = %s .\n", Y
.str());
2133 if (V
- Y
== V
+ V0
)
2134 printf ("Seems O.K.\n");
2137 printf ("finds a ");
2138 BadCond (Flaw
, "-(-Y) differs from Y.\n");
2142 BadCond (Serious
, "");
2143 printf ("overflow past %s\n\tshrinks to %s .\n", Y
.str(), Z
.str());
2147 Y
= V
* (HInvrse
* U2
- HInvrse
);
2148 Z
= Y
+ ((One
- HInvrse
) * U2
) * V
;
2158 V
= Y
* (HInvrse
* U2
- HInvrse
);
2159 V
= V
+ ((One
- HInvrse
) * U2
) * Y
;
2161 printf ("Overflow threshold is V = %s .\n", V
.str());
2163 printf ("Overflow saturates at V0 = %s .\n", V0
.str());
2165 printf ("There is no saturation value because "
2166 "the system traps on overflow.\n");
2168 printf ("No Overflow should be signaled for V * 1 = %s\n", V9
.str());
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 /*=============================================*/
2175 /*=============================================*/
2176 if (!(-V
< V
&& -V0
< V0
&& -UfThold
< V
&& UfThold
< V
))
2178 BadCond (Failure
, "Comparisons involving ");
2179 printf ("+-%s, +-%s\nand +-%s are confused by Overflow.",
2180 V
.str(), V0
.str(), UfThold
.str());
2182 /*=============================================*/
2184 /*=============================================*/
2186 for (Indx
= 1; Indx
<= 3; ++Indx
)
2204 if (Y
/ (One
- Radix
* E9
) < Z
|| Y
> (One
+ Radix
* E9
) * Z
)
2205 { /* dgh: + E9 --> * E9 */
2207 BadCond (Serious
, "");
2209 BadCond (Defect
, "");
2210 printf ("Comparison alleges that what prints as Z = %s\n",
2212 printf (" is too far from sqrt(Z) ^ 2 = %s .\n", Y
.str());
2216 /*=============================================*/
2218 /*=============================================*/
2219 for (Indx
= 1; Indx
<= 2; ++Indx
)
2226 X
= (One
- Radix
* E9
) * V9
;
2228 if (((V9
< (One
- Two
* Radix
* E9
) * Z
) || (V9
> Z
)))
2232 BadCond (Serious
, "");
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());
2239 /*=============================================*/
2241 /*=============================================*/
2245 if (X
* Y
< One
|| X
> Y
)
2247 if (X
* Y
< U1
|| X
> Y
/ U1
)
2248 BadCond (Defect
, "Badly");
2252 printf (" unbalanced range; UfThold * V = %s\n\t%s\n",
2253 X
.str(), "is too far from 1.\n");
2255 /*=============================================*/
2257 /*=============================================*/
2258 for (Indx
= 1; Indx
<= 5; ++Indx
)
2276 if (setjmp (ovfl_buf
))
2277 printf (" X / X traps when X = %s\n", X
.str());
2280 V9
= (Y
/ X
- Half
) - Half
;
2283 if (V9
== -U1
&& Indx
< 5)
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());
2291 /*=============================================*/
2293 /*=============================================*/
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 /*=============================================*/
2305 /*=============================================*/
2309 static const char *msg
[] = {
2310 "FAILUREs encountered =",
2311 "SERIOUS DEFECTs discovered =",
2312 "DEFECTs discovered =",
2313 "FLAWs discovered ="
2316 for (i
= 0; i
< 4; i
++)
2318 printf ("The number of %-29s %d.\n", msg
[i
], ErrCnt
[i
]);
2321 if ((ErrCnt
[Failure
] + ErrCnt
[Serious
] + ErrCnt
[Defect
] + ErrCnt
[Flaw
]) > 0)
2323 if ((ErrCnt
[Failure
] + ErrCnt
[Serious
] + ErrCnt
[Defect
] == 0)
2324 && (ErrCnt
[Flaw
] > 0))
2326 printf ("The arithmetic diagnosed seems ");
2327 printf ("Satisfactory though flawed.\n");
2329 if ((ErrCnt
[Failure
] + ErrCnt
[Serious
] == 0) && (ErrCnt
[Defect
] > 0))
2331 printf ("The arithmetic diagnosed may be Acceptable\n");
2332 printf ("despite inconvenient Defects.\n");
2334 if ((ErrCnt
[Failure
] + ErrCnt
[Serious
]) > 0)
2336 printf ("The arithmetic diagnosed has ");
2337 printf ("unacceptable Serious Defects.\n");
2339 if (ErrCnt
[Failure
] > 0)
2341 printf ("Potentially fatal FAILURE may have spoiled this");
2342 printf (" program's subsequent diagnoses.\n");
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");
2353 if (StickyBit
>= One
&&
2354 (Radix
- Two
) * (Radix
- Nine
- One
) == Zero
)
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
))
2368 printf (",\nexcept for possibly Double Rounding");
2369 printf (" during Gradual Underflow.\n");
2372 printf ("The arithmetic diagnosed appears to be Excellent!\n");
2375 printf ("END OF TEST.\n");
2379 template<typename FLOAT
>
2381 Paranoia
<FLOAT
>::Sign (FLOAT X
)
2383 return X
>= FLOAT (long (0)) ? 1 : -1;
2386 template<typename FLOAT
>
2388 Paranoia
<FLOAT
>::Pause ()
2392 fputs ("Press return...", stdout
);
2396 printf ("\nDiagnosis resumes after milestone Number %d", Milestone
);
2397 printf (" Page: %d\n\n", PageNo
);
2402 template<typename FLOAT
>
2404 Paranoia
<FLOAT
>::TstCond (int K
, int Valid
, const char *T
)
2413 template<typename FLOAT
>
2415 Paranoia
<FLOAT
>::BadCond (int K
, const char *T
)
2417 static const char *msg
[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
2419 ErrCnt
[K
] = ErrCnt
[K
] + 1;
2420 printf ("%s: %s", msg
[K
], T
);
2424 X = (Random1 + Random9)^5
2425 Random1 = X - FLOOR(X) + 0.000005 * X;
2426 and returns the new value of Random1. */
2428 template<typename FLOAT
>
2430 Paranoia
<FLOAT
>::Random ()
2434 X
= Random1
+ Random9
;
2439 Random1
= Y
+ X
* FLOAT ("0.000005");
2443 template<typename FLOAT
>
2445 Paranoia
<FLOAT
>::SqXMinX (int ErrKind
)
2451 SqEr
= ((SQRT (X
* X
) - XB
) - XA
) / OneUlp
;
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");
2466 template<typename FLOAT
>
2468 Paranoia
<FLOAT
>::NewD ()
2471 X
= FLOOR (Half
- X
/ Radix
) * Radix
+ X
;
2472 Q
= (Q
- X
* Z
) / Radix
+ X
* X
* (D
/ Radix
);
2473 Z
= Z
- Two
* X
* D
;
2482 template<typename FLOAT
>
2484 Paranoia
<FLOAT
>::SR3750 ()
2486 if (!((X
- Radix
< Z2
- Radix
) || (X
- Z2
> W
- Z2
)))
2490 Y2
= (X2
- Z2
) - (Y
- Z2
);
2491 X2
= X8
/ (Y
- Half
);
2492 X2
= X2
- Half
* X2
* X2
;
2493 SqEr
= (Y2
+ Half
) + (Half
- X2
);
2502 template<typename FLOAT
>
2504 Paranoia
<FLOAT
>::IsYeqX ()
2510 if (Z
== Zero
&& Q
<= Zero
)
2511 printf ("WARNING: computing\n");
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());
2519 N
= N
+ 1; /* ... count discrepancies. */
2523 template<typename FLOAT
>
2525 Paranoia
<FLOAT
>::PrintIfNPositive ()
2528 printf ("Similar discrepancies have occurred %d times.\n", N
);
2531 template<typename FLOAT
>
2533 Paranoia
<FLOAT
>::TstPtUf ()
2538 printf ("Since comparison denies Z = 0, evaluating ");
2539 printf ("(Z + Z) / Z should be safe.\n");
2540 if (setjmp (ovfl_buf
))
2543 printf ("What the machine gets for (Z + Z) / Z is %s .\n", Q9
.str());
2544 if (FABS (Q9
- Two
) < Radix
* U2
)
2546 printf ("This is O.K., provided Over/Underflow");
2547 printf (" has NOT just been signaled.\n");
2551 if ((Q9
< One
) || (Q9
> Two
))
2555 ErrCnt
[Serious
] = ErrCnt
[Serious
] + 1;
2556 printf ("This is a VERY SERIOUS DEFECT!\n");
2561 ErrCnt
[Defect
] = ErrCnt
[Defect
] + 1;
2562 printf ("This is a DEFECT!\n");
2570 if ((Z
== Random1
) && (Z
== Random2
) && (Z
== V9
))
2578 BadCond (Defect
, "What prints as Z = ");
2579 printf ("%s\n\tcompares different from ", Z
.str());
2581 printf ("Z * 1 = %s ", Random1
.str());
2582 if (!((Z
== Random2
) || (Random2
== Random1
)))
2583 printf ("1 * Z == %s\n", Random2
.str());
2585 printf ("Z / 1 = %s\n", V9
.str());
2586 if (Random2
!= Random1
)
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());
2598 template<typename FLOAT
>
2600 Paranoia
<FLOAT
>::notify (const char *s
)
2602 printf ("%s test appears to be inconsistent...\n", s
);
2603 printf (" PLEASE NOTIFY KARPINKSI!\n");
2606 /* ====================================================================== */
2608 int main(int ac
, char **av
)
2611 switch (getopt (ac
, av
, "pvg:fdl"))
2623 static const struct {
2625 const struct real_format
*fmt
;
2627 #define F(x) { #x, &x##_format }
2630 F(ieee_extended_motorola
),
2631 F(ieee_extended_intel_96
),
2632 F(ieee_extended_intel_128
),
2645 int i
, n
= sizeof (fmts
)/sizeof(*fmts
);
2647 for (i
= 0; i
< n
; ++i
)
2648 if (strcmp (fmts
[i
].name
, optarg
) == 0)
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
);
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
)]
2665 Paranoia
<real_c_float
>().main();
2670 Paranoia
< native_float
<float> >().main();
2673 Paranoia
< native_float
<double> >().main();
2676 #ifndef NO_LONG_DOUBLE
2677 Paranoia
< native_float
<long double> >().main();
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");
2691 /* GCC stuff referenced by real.o. */
2699 int target_flags
= 0;