Mercurial > hg > CbC > CbC_gcc

comparison contrib/paranoia.cc @ 0:a06113de4d67

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