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