|
0
|
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 }
|
