1 | /* |
---|
2 | * $Id$ |
---|
3 | * |
---|
4 | * A C version of Kahan's Floating Point Test "Paranoia" |
---|
5 | * |
---|
6 | * Thos Sumner, UCSF, Feb. 1985 |
---|
7 | * David Gay, BTL, Jan. 1986 |
---|
8 | * |
---|
9 | * This is a rewrite from the Pascal version by |
---|
10 | * |
---|
11 | * B. A. Wichmann, 18 Jan. 1985 |
---|
12 | * |
---|
13 | * (and does NOT exhibit good C programming style). |
---|
14 | * |
---|
15 | * Sun May 16 18:21:51 MDT 1993 Jeffrey Wheat (cassidy@cygnus.com) |
---|
16 | * Removed KR_headers defines, removed ANSI prototyping |
---|
17 | * Cleaned up and reformated code. Added special CYGNUS |
---|
18 | * "verbose" mode type messages (-DCYGNUS). |
---|
19 | * Note: This code is VERY NASTY. |
---|
20 | * |
---|
21 | * Adjusted to use Standard C headers 19 Jan. 1992 (dmg); |
---|
22 | * compile with -DKR_headers or insert |
---|
23 | * #define KR_headers |
---|
24 | * at the beginning if you have an old-style C compiler. |
---|
25 | * |
---|
26 | * (C) Apr 19 1983 in BASIC version by: |
---|
27 | * Professor W. M. Kahan, |
---|
28 | * 567 Evans Hall |
---|
29 | * Electrical Engineering & Computer Science Dept. |
---|
30 | * University of California |
---|
31 | * Berkeley, California 94720 |
---|
32 | * USA |
---|
33 | * |
---|
34 | * converted to Pascal by: |
---|
35 | * B. A. Wichmann |
---|
36 | * National Physical Laboratory |
---|
37 | * Teddington Middx |
---|
38 | * TW11 OLW |
---|
39 | * UK |
---|
40 | * |
---|
41 | * converted to C by: |
---|
42 | * |
---|
43 | * David M. Gay and Thos Sumner |
---|
44 | * AT&T Bell Labs Computer Center, Rm. U-76 |
---|
45 | * 600 Mountain Avenue University of California |
---|
46 | * Murray Hill, NJ 07974 San Francisco, CA 94143 |
---|
47 | * USA USA |
---|
48 | * |
---|
49 | * with simultaneous corrections to the Pascal source (reflected |
---|
50 | * in the Pascal source available over netlib). |
---|
51 | * [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.] |
---|
52 | * |
---|
53 | * Reports of results on various systems from all the versions |
---|
54 | * of Paranoia are being collected by Richard Karpinski at the |
---|
55 | * same address as Thos Sumner. This includes sample outputs, |
---|
56 | * bug reports, and criticisms. |
---|
57 | * |
---|
58 | * You may copy this program freely if you acknowledge its source. |
---|
59 | * Comments on the Pascal version to NPL, please. |
---|
60 | * |
---|
61 | * |
---|
62 | * The C version catches signals from floating-point exceptions. |
---|
63 | * If signal(SIGFPE,...) is unavailable in your environment, you may |
---|
64 | * #define NOSIGNAL to comment out the invocations of signal. |
---|
65 | * |
---|
66 | * This source file is too big for some C compilers, but may be split |
---|
67 | * into pieces. Comments containing "SPLIT" suggest convenient places |
---|
68 | * for this splitting. At the end of these comments is an "ed script" |
---|
69 | * (for the UNIX(tm) editor ed) that will do this splitting. |
---|
70 | * |
---|
71 | * By #defining SINGLE_PRECISION when you compile this source, you may |
---|
72 | * obtain a single-precision C version of Paranoia. |
---|
73 | * |
---|
74 | * The following is from the introductory commentary from Wichmann's work: |
---|
75 | * |
---|
76 | * The BASIC program of Kahan is written in Microsoft BASIC using many |
---|
77 | * facilities which have no exact analogy in Pascal. The Pascal |
---|
78 | * version below cannot therefore be exactly the same. Rather than be |
---|
79 | * a minimal transcription of the BASIC program, the Pascal coding |
---|
80 | * follows the conventional style of block-structured languages. Hence |
---|
81 | * the Pascal version could be useful in producing versions in other |
---|
82 | * structured languages. |
---|
83 | * |
---|
84 | * Rather than use identifiers of minimal length (which therefore have |
---|
85 | * little mnemonic significance), the Pascal version uses meaningful |
---|
86 | * identifiers as follows [Note: A few changes have been made for C]: |
---|
87 | * |
---|
88 | * |
---|
89 | * BASIC C BASIC C BASIC C |
---|
90 | * |
---|
91 | * A J S StickyBit |
---|
92 | * A1 AInverse J0 NoErrors T |
---|
93 | * B Radix [Failure] T0 Underflow |
---|
94 | * B1 BInverse J1 NoErrors T2 ThirtyTwo |
---|
95 | * B2 RadixD2 [SeriousDefect] T5 OneAndHalf |
---|
96 | * B9 BMinusU2 J2 NoErrors T7 TwentySeven |
---|
97 | * C [Defect] T8 TwoForty |
---|
98 | * C1 CInverse J3 NoErrors U OneUlp |
---|
99 | * D [Flaw] U0 UnderflowThreshold |
---|
100 | * D4 FourD K PageNo U1 |
---|
101 | * E0 L Milestone U2 |
---|
102 | * E1 M V |
---|
103 | * E2 Exp2 N V0 |
---|
104 | * E3 N1 V8 |
---|
105 | * E5 MinSqEr O Zero V9 |
---|
106 | * E6 SqEr O1 One W |
---|
107 | * E7 MaxSqEr O2 Two X |
---|
108 | * E8 O3 Three X1 |
---|
109 | * E9 O4 Four X8 |
---|
110 | * F1 MinusOne O5 Five X9 Random1 |
---|
111 | * F2 Half O8 Eight Y |
---|
112 | * F3 Third O9 Nine Y1 |
---|
113 | * F6 P Precision Y2 |
---|
114 | * F9 Q Y9 Random2 |
---|
115 | * G1 GMult Q8 Z |
---|
116 | * G2 GDiv Q9 Z0 PseudoZero |
---|
117 | * G3 GAddSub R Z1 |
---|
118 | * H R1 RMult Z2 |
---|
119 | * H1 HInverse R2 RDiv Z9 |
---|
120 | * I R3 RAddSub |
---|
121 | * IO NoTrials R4 RSqrt |
---|
122 | * I3 IEEE R9 Random9 |
---|
123 | * |
---|
124 | * SqRWrng |
---|
125 | * |
---|
126 | * All the variables in BASIC are true variables and in consequence, |
---|
127 | * the program is more difficult to follow since the "constants" must |
---|
128 | * be determined (the glossary is very helpful). The Pascal version |
---|
129 | * uses Real constants, but checks are added to ensure that the values |
---|
130 | * are correctly converted by the compiler. |
---|
131 | * |
---|
132 | * The major textual change to the Pascal version apart from the |
---|
133 | * identifiersis that named procedures are used, inserting parameters |
---|
134 | * wherehelpful. New procedures are also introduced. The |
---|
135 | * correspondence is as follows: |
---|
136 | * |
---|
137 | * |
---|
138 | * BASIC Pascal |
---|
139 | * lines |
---|
140 | * |
---|
141 | * 90- 140 Pause |
---|
142 | * 170- 250 Instructions |
---|
143 | * 380- 460 Heading |
---|
144 | * 480- 670 Characteristics |
---|
145 | * 690- 870 History |
---|
146 | * 2940-2950 Random |
---|
147 | * 3710-3740 NewD |
---|
148 | * 4040-4080 DoesYequalX |
---|
149 | * 4090-4110 PrintIfNPositive |
---|
150 | * 4640-4850 TestPartialUnderflow |
---|
151 | * |
---|
152 | */ |
---|
153 | |
---|
154 | #include <stdio.h> |
---|
155 | #include <string.h> |
---|
156 | |
---|
157 | #if defined(solaris2) |
---|
158 | #include <math.h> |
---|
159 | #endif |
---|
160 | |
---|
161 | /* |
---|
162 | * To compile this on host using only libm from newlib (and using host libc) |
---|
163 | * do: |
---|
164 | * gcc -g -DNEED_REENT -DCYGNUS paranoia.c .../newlib-1.6/newlib/libm.a |
---|
165 | */ |
---|
166 | |
---|
167 | #ifdef NEED_REENT |
---|
168 | #include <reent.h> |
---|
169 | struct _reent libm_reent = _REENT_INIT(libm_reent); |
---|
170 | struct _reent *_impure_ptr = &libm_reent; |
---|
171 | #endif |
---|
172 | |
---|
173 | #ifndef NOSIGNAL |
---|
174 | #include <signal.h> |
---|
175 | #include <setjmp.h> |
---|
176 | #else /* NOSIGNAL */ |
---|
177 | #define longjmp(e,v) |
---|
178 | #define setjmp(e) 0 |
---|
179 | #define jmp_buf int |
---|
180 | #endif /* NOSIGNAL */ |
---|
181 | |
---|
182 | #ifdef SINGLE_PRECISION |
---|
183 | #define FLOAT float |
---|
184 | #define FABS(x) (float)fabs((double)(x)) |
---|
185 | #define FLOOR(x) (float)floor((double)(x)) |
---|
186 | #define LOG(x) (float)log((double)(x)) |
---|
187 | #define POW(x,y) (float)pow((double)(x),(double)(y)) |
---|
188 | #define SQRT(x) (float)sqrt((double)(x)) |
---|
189 | #else /* !SINGLE_PRECISION */ |
---|
190 | #define FLOAT double |
---|
191 | #define FABS(x) fabs(x) |
---|
192 | #define FLOOR(x) floor(x) |
---|
193 | #define LOG(x) log(x) |
---|
194 | #define POW(x,y) pow(x,y) |
---|
195 | #define SQRT(x) sqrt(x) |
---|
196 | #endif /* SINGLE_PRECISION */ |
---|
197 | |
---|
198 | jmp_buf ovfl_buf; |
---|
199 | extern double fabs (), floor (), log (), pow (), sqrt (); |
---|
200 | extern void exit (); |
---|
201 | typedef void (*Sig_type) (); |
---|
202 | FLOAT Sign (), Random (); |
---|
203 | extern void BadCond (); |
---|
204 | extern void SqXMinX (); |
---|
205 | extern void TstCond (); |
---|
206 | extern void notify (); |
---|
207 | extern int read (); |
---|
208 | extern void Characteristics (); |
---|
209 | extern void Heading (); |
---|
210 | extern void History (); |
---|
211 | extern void Instructions (); |
---|
212 | extern void IsYeqX (); |
---|
213 | extern void NewD (); |
---|
214 | extern void Pause (); |
---|
215 | extern void PrintIfNPositive (); |
---|
216 | extern void SR3750 (); |
---|
217 | extern void SR3980 (); |
---|
218 | extern void TstPtUf (); |
---|
219 | |
---|
220 | Sig_type sigsave; |
---|
221 | |
---|
222 | #define KEYBOARD 0 |
---|
223 | |
---|
224 | FLOAT Radix, BInvrse, RadixD2, BMinusU2; |
---|
225 | |
---|
226 | /*Small floating point constants.*/ |
---|
227 | FLOAT Zero = 0.0; |
---|
228 | FLOAT Half = 0.5; |
---|
229 | FLOAT One = 1.0; |
---|
230 | FLOAT Two = 2.0; |
---|
231 | FLOAT Three = 3.0; |
---|
232 | FLOAT Four = 4.0; |
---|
233 | FLOAT Five = 5.0; |
---|
234 | FLOAT Eight = 8.0; |
---|
235 | FLOAT Nine = 9.0; |
---|
236 | FLOAT TwentySeven = 27.0; |
---|
237 | FLOAT ThirtyTwo = 32.0; |
---|
238 | FLOAT TwoForty = 240.0; |
---|
239 | FLOAT MinusOne = -1.0; |
---|
240 | FLOAT OneAndHalf = 1.5; |
---|
241 | |
---|
242 | /*Integer constants*/ |
---|
243 | int NoTrials = 20; /*Number of tests for commutativity. */ |
---|
244 | #define False 0 |
---|
245 | #define True 1 |
---|
246 | |
---|
247 | /* |
---|
248 | * Definitions for declared types |
---|
249 | * Guard == (Yes, No); |
---|
250 | * Rounding == (Chopped, Rounded, Other); |
---|
251 | * Message == packed array [1..40] of char; |
---|
252 | * Class == (Flaw, Defect, Serious, Failure); |
---|
253 | */ |
---|
254 | #define Yes 1 |
---|
255 | #define No 0 |
---|
256 | #define Chopped 2 |
---|
257 | #define Rounded 1 |
---|
258 | #define Other 0 |
---|
259 | #define Flaw 3 |
---|
260 | #define Defect 2 |
---|
261 | #define Serious 1 |
---|
262 | #define Failure 0 |
---|
263 | typedef int Guard, Rounding, Class; |
---|
264 | typedef char Message; |
---|
265 | |
---|
266 | /* Declarations of Variables */ |
---|
267 | int Indx; |
---|
268 | char ch[8]; |
---|
269 | FLOAT AInvrse, A1; |
---|
270 | FLOAT C, CInvrse; |
---|
271 | FLOAT D, FourD; |
---|
272 | FLOAT E0, E1, Exp2, E3, MinSqEr; |
---|
273 | FLOAT SqEr, MaxSqEr, E9; |
---|
274 | FLOAT Third; |
---|
275 | FLOAT F6, F9; |
---|
276 | FLOAT HVar, HInvrse; |
---|
277 | int I; |
---|
278 | FLOAT StickyBit, J; |
---|
279 | FLOAT MyZero; |
---|
280 | FLOAT Precision; |
---|
281 | FLOAT Q, Q9; |
---|
282 | FLOAT R, Random9; |
---|
283 | FLOAT T, Underflow, S; |
---|
284 | FLOAT OneUlp, UfThold, U1, U2; |
---|
285 | FLOAT V, V0, V9; |
---|
286 | FLOAT WVar; |
---|
287 | FLOAT X, X1, X2, X8, Random1; |
---|
288 | FLOAT Y, Y1, Y2, Random2; |
---|
289 | FLOAT Z, PseudoZero, Z1, Z2, Z9; |
---|
290 | int ErrCnt[4]; |
---|
291 | int fpecount; |
---|
292 | int Milestone; |
---|
293 | int PageNo; |
---|
294 | int M, N, N1; |
---|
295 | Guard GMult, GDiv, GAddSub; |
---|
296 | Rounding RMult, RDiv, RAddSub, RSqrt; |
---|
297 | int Break, Done, NotMonot, Monot, Anomaly, IEEE, SqRWrng, UfNGrad; |
---|
298 | |
---|
299 | /* Computed constants. |
---|
300 | * U1 gap below 1.0, i.e, 1.0 - U1 is next number below 1.0 |
---|
301 | * U2 gap above 1.0, i.e, 1.0 + U2 is next number above 1.0 |
---|
302 | */ |
---|
303 | |
---|
304 | int batchmode; /* global batchmode test */ |
---|
305 | |
---|
306 | /* program name and version variables and macro */ |
---|
307 | char *temp; |
---|
308 | char *program_name; |
---|
309 | char *program_vers; |
---|
310 | |
---|
311 | #ifndef VERSION |
---|
312 | #define VERSION "1.1 [cygnus]" |
---|
313 | #endif /* VERSION */ |
---|
314 | |
---|
315 | #define basename(cp) ((temp=(char *)strrchr((cp), '/')) ? temp+1 : (cp)) |
---|
316 | |
---|
317 | #ifndef BATCHMODE |
---|
318 | # ifdef CYGNUS |
---|
319 | # define BATCHMODE |
---|
320 | # endif |
---|
321 | #endif |
---|
322 | |
---|
323 | /* floating point exception receiver */ |
---|
324 | void |
---|
325 | _sigfpe (x) |
---|
326 | int x; |
---|
327 | { |
---|
328 | fpecount++; |
---|
329 | printf ("\n* * * FLOATING-POINT ERROR %d * * *\n", x); |
---|
330 | fflush (stdout); |
---|
331 | if (sigsave) { |
---|
332 | #ifndef NOSIGNAL |
---|
333 | signal (SIGFPE, sigsave); |
---|
334 | #endif /* NOSIGNAL */ |
---|
335 | sigsave = 0; |
---|
336 | longjmp (ovfl_buf, 1); |
---|
337 | } |
---|
338 | exit (1); |
---|
339 | } |
---|
340 | |
---|
341 | #ifdef NOMAIN |
---|
342 | #define main paranoia |
---|
343 | #endif |
---|
344 | |
---|
345 | int |
---|
346 | main (argc, argv) |
---|
347 | int argc; |
---|
348 | char **argv; |
---|
349 | { |
---|
350 | /* First two assignments use integer right-hand sides. */ |
---|
351 | Zero = 0; |
---|
352 | One = 1; |
---|
353 | Two = One + One; |
---|
354 | Three = Two + One; |
---|
355 | Four = Three + One; |
---|
356 | Five = Four + One; |
---|
357 | Eight = Four + Four; |
---|
358 | Nine = Three * Three; |
---|
359 | TwentySeven = Nine * Three; |
---|
360 | ThirtyTwo = Four * Eight; |
---|
361 | TwoForty = Four * Five * Three * Four; |
---|
362 | MinusOne = -One; |
---|
363 | Half = One / Two; |
---|
364 | OneAndHalf = One + Half; |
---|
365 | ErrCnt[Failure] = 0; |
---|
366 | ErrCnt[Serious] = 0; |
---|
367 | ErrCnt[Defect] = 0; |
---|
368 | ErrCnt[Flaw] = 0; |
---|
369 | PageNo = 1; |
---|
370 | |
---|
371 | #ifdef BATCHMODE |
---|
372 | batchmode = 1; /* run test in batchmode? */ |
---|
373 | #else /* !BATCHMODE */ |
---|
374 | batchmode = 0; /* run test interactively */ |
---|
375 | #endif /* BATCHMODE */ |
---|
376 | |
---|
377 | program_name = basename (argv[0]); |
---|
378 | program_vers = VERSION; |
---|
379 | |
---|
380 | printf ("%s version %s\n", program_name, program_vers); |
---|
381 | |
---|
382 | /*=============================================*/ |
---|
383 | Milestone = 0; |
---|
384 | /*=============================================*/ |
---|
385 | #ifndef NOSIGNAL |
---|
386 | signal (SIGFPE, _sigfpe); |
---|
387 | #endif |
---|
388 | |
---|
389 | if (!batchmode) { |
---|
390 | Instructions (); |
---|
391 | Pause (); |
---|
392 | Heading (); |
---|
393 | Instructions (); |
---|
394 | Pause (); |
---|
395 | Heading (); |
---|
396 | Pause (); |
---|
397 | Characteristics (); |
---|
398 | Pause (); |
---|
399 | History (); |
---|
400 | Pause (); |
---|
401 | } |
---|
402 | |
---|
403 | /*=============================================*/ |
---|
404 | Milestone = 7; |
---|
405 | /*=============================================*/ |
---|
406 | printf ("Program is now RUNNING tests on small integers:\n"); |
---|
407 | TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero) |
---|
408 | && (One > Zero) && (One + One == Two), |
---|
409 | "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2"); |
---|
410 | Z = -Zero; |
---|
411 | if (Z != 0.0) { |
---|
412 | ErrCnt[Failure] = ErrCnt[Failure] + 1; |
---|
413 | printf ("Comparison alleges that -0.0 is Non-zero!\n"); |
---|
414 | U1 = 0.001; |
---|
415 | Radix = 1; |
---|
416 | TstPtUf (); |
---|
417 | } |
---|
418 | TstCond (Failure, (Three == Two + One) && (Four == Three + One) |
---|
419 | && (Four + Two * (-Two) == Zero) |
---|
420 | && (Four - Three - One == Zero), |
---|
421 | "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0"); |
---|
422 | TstCond (Failure, (MinusOne == (0 - One)) |
---|
423 | && (MinusOne + One == Zero) && (One + MinusOne == Zero) |
---|
424 | && (MinusOne + FABS (One) == Zero) |
---|
425 | && (MinusOne + MinusOne * MinusOne == Zero), |
---|
426 | "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0"); |
---|
427 | TstCond (Failure, Half + MinusOne + Half == Zero, |
---|
428 | "1/2 + (-1) + 1/2 != 0"); |
---|
429 | /*=============================================*/ |
---|
430 | Milestone = 10; |
---|
431 | /*=============================================*/ |
---|
432 | TstCond (Failure, (Nine == Three * Three) |
---|
433 | && (TwentySeven == Nine * Three) && (Eight == Four + Four) |
---|
434 | && (ThirtyTwo == Eight * Four) |
---|
435 | && (ThirtyTwo - TwentySeven - Four - One == Zero), |
---|
436 | "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0"); |
---|
437 | TstCond (Failure, (Five == Four + One) && |
---|
438 | (TwoForty == Four * Five * Three * Four) |
---|
439 | && (TwoForty / Three - Four * Four * Five == Zero) |
---|
440 | && (TwoForty / Four - Five * Three * Four == Zero) |
---|
441 | && (TwoForty / Five - Four * Three * Four == Zero), |
---|
442 | "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48"); |
---|
443 | if (ErrCnt[Failure] == 0) { |
---|
444 | printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n"); |
---|
445 | printf ("\n"); |
---|
446 | } |
---|
447 | printf ("Searching for Radix and Precision.\n"); |
---|
448 | WVar = One; |
---|
449 | do { |
---|
450 | WVar = WVar + WVar; |
---|
451 | Y = WVar + One; |
---|
452 | Z = Y - WVar; |
---|
453 | Y = Z - One; |
---|
454 | } |
---|
455 | while (MinusOne + FABS (Y) < Zero); |
---|
456 | /*.. now WVar is just big enough that |((WVar+1)-WVar)-1| >= 1 ...*/ |
---|
457 | Precision = Zero; |
---|
458 | Y = One; |
---|
459 | do { |
---|
460 | Radix = WVar + Y; |
---|
461 | Y = Y + Y; |
---|
462 | Radix = Radix - WVar; |
---|
463 | } |
---|
464 | while (Radix == Zero); |
---|
465 | if (Radix < Two) |
---|
466 | Radix = One; |
---|
467 | printf ("Radix = %f .\n", Radix); |
---|
468 | if (Radix != 1) { |
---|
469 | WVar = One; |
---|
470 | do { |
---|
471 | Precision = Precision + One; |
---|
472 | WVar = WVar * Radix; |
---|
473 | Y = WVar + One; |
---|
474 | } |
---|
475 | while ((Y - WVar) == One); |
---|
476 | } |
---|
477 | /*... now WVar == Radix^Precision is barely too big to satisfy (WVar+1)-WVar == 1 |
---|
478 | ...*/ |
---|
479 | U1 = One / WVar; |
---|
480 | U2 = Radix * U1; |
---|
481 | printf ("Closest relative separation found is U1 = %.7e .\n\n", U1); |
---|
482 | printf ("Recalculating radix and precision\n "); |
---|
483 | |
---|
484 | /*save old values*/ |
---|
485 | E0 = Radix; |
---|
486 | E1 = U1; |
---|
487 | E9 = U2; |
---|
488 | E3 = Precision; |
---|
489 | |
---|
490 | X = Four / Three; |
---|
491 | Third = X - One; |
---|
492 | F6 = Half - Third; |
---|
493 | X = F6 + F6; |
---|
494 | X = FABS (X - Third); |
---|
495 | if (X < U2) |
---|
496 | X = U2; |
---|
497 | |
---|
498 | /*... now X = (unknown no.) ulps of 1+...*/ |
---|
499 | do { |
---|
500 | U2 = X; |
---|
501 | Y = Half * U2 + ThirtyTwo * U2 * U2; |
---|
502 | Y = One + Y; |
---|
503 | X = Y - One; |
---|
504 | } |
---|
505 | while (!((U2 <= X) || (X <= Zero))); |
---|
506 | |
---|
507 | /*... now U2 == 1 ulp of 1 + ... */ |
---|
508 | X = Two / Three; |
---|
509 | F6 = X - Half; |
---|
510 | Third = F6 + F6; |
---|
511 | X = Third - Half; |
---|
512 | X = FABS (X + F6); |
---|
513 | if (X < U1) |
---|
514 | X = U1; |
---|
515 | |
---|
516 | /*... now X == (unknown no.) ulps of 1 -... */ |
---|
517 | do { |
---|
518 | U1 = X; |
---|
519 | Y = Half * U1 + ThirtyTwo * U1 * U1; |
---|
520 | Y = Half - Y; |
---|
521 | X = Half + Y; |
---|
522 | Y = Half - X; |
---|
523 | X = Half + Y; |
---|
524 | } |
---|
525 | while (!((U1 <= X) || (X <= Zero))); |
---|
526 | /*... now U1 == 1 ulp of 1 - ... */ |
---|
527 | if (U1 == E1) |
---|
528 | printf ("confirms closest relative separation U1 .\n"); |
---|
529 | else |
---|
530 | printf ("gets better closest relative separation U1 = %.7e .\n", U1); |
---|
531 | WVar = One / U1; |
---|
532 | F9 = (Half - U1) + Half; |
---|
533 | Radix = FLOOR (0.01 + U2 / U1); |
---|
534 | if (Radix == E0) |
---|
535 | printf ("Radix confirmed.\n"); |
---|
536 | else |
---|
537 | printf ("MYSTERY: recalculated Radix = %.7e .\n", Radix); |
---|
538 | TstCond (Defect, Radix <= Eight + Eight, |
---|
539 | "Radix is too big: roundoff problems"); |
---|
540 | TstCond (Flaw, (Radix == Two) || (Radix == 10) |
---|
541 | || (Radix == One), "Radix is not as good as 2 or 10"); |
---|
542 | /*=============================================*/ |
---|
543 | Milestone = 20; |
---|
544 | /*=============================================*/ |
---|
545 | TstCond (Failure, F9 - Half < Half, |
---|
546 | "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?"); |
---|
547 | X = F9; |
---|
548 | I = 1; |
---|
549 | Y = X - Half; |
---|
550 | Z = Y - Half; |
---|
551 | TstCond (Failure, (X != One) |
---|
552 | || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0"); |
---|
553 | X = One + U2; |
---|
554 | I = 0; |
---|
555 | /*=============================================*/ |
---|
556 | Milestone = 25; |
---|
557 | /*=============================================*/ |
---|
558 | /*... BMinusU2 = nextafter(Radix, 0) */ |
---|
559 | BMinusU2 = Radix - One; |
---|
560 | BMinusU2 = (BMinusU2 - U2) + One; |
---|
561 | /* Purify Integers */ |
---|
562 | if (Radix != One) { |
---|
563 | X = -TwoForty * LOG (U1) / LOG (Radix); |
---|
564 | Y = FLOOR (Half + X); |
---|
565 | if (FABS (X - Y) * Four < One) |
---|
566 | X = Y; |
---|
567 | Precision = X / TwoForty; |
---|
568 | Y = FLOOR (Half + Precision); |
---|
569 | if (FABS (Precision - Y) * TwoForty < Half) |
---|
570 | Precision = Y; |
---|
571 | } |
---|
572 | if ((Precision != FLOOR (Precision)) || (Radix == One)) { |
---|
573 | printf ("Precision cannot be characterized by an Integer number\n"); |
---|
574 | printf ("of significant digits but, by itself, this is a minor flaw.\n"); |
---|
575 | } |
---|
576 | if (Radix == One) |
---|
577 | printf ("logarithmic encoding has precision characterized solely by U1.\n"); |
---|
578 | else |
---|
579 | printf ("The number of significant digits of the Radix is %f .\n", |
---|
580 | Precision); |
---|
581 | TstCond (Serious, U2 * Nine * Nine * TwoForty < One, |
---|
582 | "Precision worse than 5 decimal figures "); |
---|
583 | /*=============================================*/ |
---|
584 | Milestone = 30; |
---|
585 | /*=============================================*/ |
---|
586 | /* Test for extra-precise subepressions */ |
---|
587 | X = FABS (((Four / Three - One) - One / Four) * Three - One / Four); |
---|
588 | do { |
---|
589 | Z2 = X; |
---|
590 | X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One; |
---|
591 | } |
---|
592 | while (!((Z2 <= X) || (X <= Zero))); |
---|
593 | X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four); |
---|
594 | do { |
---|
595 | Z1 = Z; |
---|
596 | Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1)) |
---|
597 | + One / Two)) + One / Two; |
---|
598 | } |
---|
599 | while (!((Z1 <= Z) || (Z <= Zero))); |
---|
600 | do { |
---|
601 | do { |
---|
602 | Y1 = Y; |
---|
603 | Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half |
---|
604 | )) + Half; |
---|
605 | } |
---|
606 | while (!((Y1 <= Y) || (Y <= Zero))); |
---|
607 | X1 = X; |
---|
608 | X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9; |
---|
609 | } |
---|
610 | while (!((X1 <= X) || (X <= Zero))); |
---|
611 | if ((X1 != Y1) || (X1 != Z1)) { |
---|
612 | BadCond (Serious, "Disagreements among the values X1, Y1, Z1,\n"); |
---|
613 | printf ("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1); |
---|
614 | printf ("are symptoms of inconsistencies introduced\n"); |
---|
615 | printf ("by extra-precise evaluation of arithmetic subexpressions.\n"); |
---|
616 | notify ("Possibly some part of this"); |
---|
617 | if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) |
---|
618 | printf ( |
---|
619 | "That feature is not tested further by this program.\n"); |
---|
620 | } else { |
---|
621 | if ((Z1 != U1) || (Z2 != U2)) { |
---|
622 | if ((Z1 >= U1) || (Z2 >= U2)) { |
---|
623 | BadCond (Failure, ""); |
---|
624 | notify ("Precision"); |
---|
625 | printf ("\tU1 = %.7e, Z1 - U1 = %.7e\n", U1, Z1 - U1); |
---|
626 | printf ("\tU2 = %.7e, Z2 - U2 = %.7e\n", U2, Z2 - U2); |
---|
627 | } else { |
---|
628 | if ((Z1 <= Zero) || (Z2 <= Zero)) { |
---|
629 | printf ("Because of unusual Radix = %f", Radix); |
---|
630 | printf (", or exact rational arithmetic a result\n"); |
---|
631 | printf ("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2); |
---|
632 | notify ("of an\nextra-precision"); |
---|
633 | } |
---|
634 | if (Z1 != Z2 || Z1 > Zero) { |
---|
635 | X = Z1 / U1; |
---|
636 | Y = Z2 / U2; |
---|
637 | if (Y > X) |
---|
638 | X = Y; |
---|
639 | Q = -LOG (X); |
---|
640 | printf ("Some subexpressions appear to be calculated extra\n"); |
---|
641 | printf ("precisely with about %g extra B-digits, i.e.\n", |
---|
642 | (Q / LOG (Radix))); |
---|
643 | printf ("roughly %g extra significant decimals.\n", |
---|
644 | Q / LOG (10.)); |
---|
645 | } |
---|
646 | printf ("That feature is not tested further by this program.\n"); |
---|
647 | } |
---|
648 | } |
---|
649 | } |
---|
650 | Pause (); |
---|
651 | /*=============================================*/ |
---|
652 | Milestone = 35; |
---|
653 | /*=============================================*/ |
---|
654 | if (Radix >= Two) { |
---|
655 | X = WVar / (Radix * Radix); |
---|
656 | Y = X + One; |
---|
657 | Z = Y - X; |
---|
658 | T = Z + U2; |
---|
659 | X = T - Z; |
---|
660 | TstCond (Failure, X == U2, |
---|
661 | "Subtraction is not normalized X=Y,X+Z != Y+Z!"); |
---|
662 | if (X == U2) |
---|
663 | printf ( |
---|
664 | "Subtraction appears to be normalized, as it should be."); |
---|
665 | } |
---|
666 | printf ("\nChecking for guard digit in *, /, and -.\n"); |
---|
667 | Y = F9 * One; |
---|
668 | Z = One * F9; |
---|
669 | X = F9 - Half; |
---|
670 | Y = (Y - Half) - X; |
---|
671 | Z = (Z - Half) - X; |
---|
672 | X = One + U2; |
---|
673 | T = X * Radix; |
---|
674 | R = Radix * X; |
---|
675 | X = T - Radix; |
---|
676 | X = X - Radix * U2; |
---|
677 | T = R - Radix; |
---|
678 | T = T - Radix * U2; |
---|
679 | X = X * (Radix - One); |
---|
680 | T = T * (Radix - One); |
---|
681 | if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) |
---|
682 | GMult = Yes; |
---|
683 | else { |
---|
684 | GMult = No; |
---|
685 | TstCond (Serious, False, |
---|
686 | "* lacks a Guard Digit, so 1*X != X"); |
---|
687 | } |
---|
688 | Z = Radix * U2; |
---|
689 | X = One + Z; |
---|
690 | Y = FABS ((X + Z) - X * X) - U2; |
---|
691 | X = One - U2; |
---|
692 | Z = FABS ((X - U2) - X * X) - U1; |
---|
693 | TstCond (Failure, (Y <= Zero) |
---|
694 | && (Z <= Zero), "* gets too many final digits wrong.\n"); |
---|
695 | Y = One - U2; |
---|
696 | X = One + U2; |
---|
697 | Z = One / Y; |
---|
698 | Y = Z - X; |
---|
699 | X = One / Three; |
---|
700 | Z = Three / Nine; |
---|
701 | X = X - Z; |
---|
702 | T = Nine / TwentySeven; |
---|
703 | Z = Z - T; |
---|
704 | TstCond (Defect, X == Zero && Y == Zero && Z == Zero, |
---|
705 | "Division lacks a Guard Digit, so error can exceed 1 ulp\n\ |
---|
706 | or 1/3 and 3/9 and 9/27 may disagree"); |
---|
707 | Y = F9 / One; |
---|
708 | X = F9 - Half; |
---|
709 | Y = (Y - Half) - X; |
---|
710 | X = One + U2; |
---|
711 | T = X / One; |
---|
712 | X = T - X; |
---|
713 | if ((X == Zero) && (Y == Zero) && (Z == Zero)) |
---|
714 | GDiv = Yes; |
---|
715 | else { |
---|
716 | GDiv = No; |
---|
717 | TstCond (Serious, False, |
---|
718 | "Division lacks a Guard Digit, so X/1 != X"); |
---|
719 | } |
---|
720 | X = One / (One + U2); |
---|
721 | Y = X - Half - Half; |
---|
722 | TstCond (Serious, Y < Zero, |
---|
723 | "Computed value of 1/1.000..1 >= 1"); |
---|
724 | X = One - U2; |
---|
725 | Y = One + Radix * U2; |
---|
726 | Z = X * Radix; |
---|
727 | T = Y * Radix; |
---|
728 | R = Z / Radix; |
---|
729 | StickyBit = T / Radix; |
---|
730 | X = R - X; |
---|
731 | Y = StickyBit - Y; |
---|
732 | TstCond (Failure, X == Zero && Y == Zero, |
---|
733 | "* and/or / gets too many last digits wrong"); |
---|
734 | Y = One - U1; |
---|
735 | X = One - F9; |
---|
736 | Y = One - Y; |
---|
737 | T = Radix - U2; |
---|
738 | Z = Radix - BMinusU2; |
---|
739 | T = Radix - T; |
---|
740 | if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) |
---|
741 | GAddSub = Yes; |
---|
742 | else { |
---|
743 | GAddSub = No; |
---|
744 | TstCond (Serious, False, |
---|
745 | "- lacks Guard Digit, so cancellation is obscured"); |
---|
746 | } |
---|
747 | if (F9 != One && F9 - One >= Zero) { |
---|
748 | BadCond (Serious, "comparison alleges (1-U1) < 1 although\n"); |
---|
749 | printf (" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n"); |
---|
750 | printf (" such precautions against division by zero as\n"); |
---|
751 | printf (" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n"); |
---|
752 | } |
---|
753 | if (GMult == Yes && GDiv == Yes && GAddSub == Yes) |
---|
754 | printf ( |
---|
755 | " *, /, and - appear to have guard digits, as they should.\n"); |
---|
756 | /*=============================================*/ |
---|
757 | Milestone = 40; |
---|
758 | /*=============================================*/ |
---|
759 | Pause (); |
---|
760 | printf ("Checking rounding on multiply, divide and add/subtract.\n"); |
---|
761 | RMult = Other; |
---|
762 | RDiv = Other; |
---|
763 | RAddSub = Other; |
---|
764 | RadixD2 = Radix / Two; |
---|
765 | A1 = Two; |
---|
766 | Done = False; |
---|
767 | do { |
---|
768 | AInvrse = Radix; |
---|
769 | do { |
---|
770 | X = AInvrse; |
---|
771 | AInvrse = AInvrse / A1; |
---|
772 | } |
---|
773 | while (!(FLOOR (AInvrse) != AInvrse)); |
---|
774 | Done = (X == One) || (A1 > Three); |
---|
775 | if (!Done) |
---|
776 | A1 = Nine + One; |
---|
777 | } |
---|
778 | while (!(Done)); |
---|
779 | if (X == One) |
---|
780 | A1 = Radix; |
---|
781 | AInvrse = One / A1; |
---|
782 | X = A1; |
---|
783 | Y = AInvrse; |
---|
784 | Done = False; |
---|
785 | do { |
---|
786 | Z = X * Y - Half; |
---|
787 | TstCond (Failure, Z == Half, |
---|
788 | "X * (1/X) differs from 1"); |
---|
789 | Done = X == Radix; |
---|
790 | X = Radix; |
---|
791 | Y = One / X; |
---|
792 | } |
---|
793 | while (!(Done)); |
---|
794 | Y2 = One + U2; |
---|
795 | Y1 = One - U2; |
---|
796 | X = OneAndHalf - U2; |
---|
797 | Y = OneAndHalf + U2; |
---|
798 | Z = (X - U2) * Y2; |
---|
799 | T = Y * Y1; |
---|
800 | Z = Z - X; |
---|
801 | T = T - X; |
---|
802 | X = X * Y2; |
---|
803 | Y = (Y + U2) * Y1; |
---|
804 | X = X - OneAndHalf; |
---|
805 | Y = Y - OneAndHalf; |
---|
806 | if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) { |
---|
807 | X = (OneAndHalf + U2) * Y2; |
---|
808 | Y = OneAndHalf - U2 - U2; |
---|
809 | Z = OneAndHalf + U2 + U2; |
---|
810 | T = (OneAndHalf - U2) * Y1; |
---|
811 | X = X - (Z + U2); |
---|
812 | StickyBit = Y * Y1; |
---|
813 | S = Z * Y2; |
---|
814 | T = T - Y; |
---|
815 | Y = (U2 - Y) + StickyBit; |
---|
816 | Z = S - (Z + U2 + U2); |
---|
817 | StickyBit = (Y2 + U2) * Y1; |
---|
818 | Y1 = Y2 * Y1; |
---|
819 | StickyBit = StickyBit - Y2; |
---|
820 | Y1 = Y1 - Half; |
---|
821 | if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero) |
---|
822 | && (StickyBit == Zero) && (Y1 == Half)) { |
---|
823 | RMult = Rounded; |
---|
824 | printf ("Multiplication appears to round correctly.\n"); |
---|
825 | } else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero) |
---|
826 | && (T < Zero) && (StickyBit + U2 == Zero) |
---|
827 | && (Y1 < Half)) { |
---|
828 | RMult = Chopped; |
---|
829 | printf ("Multiplication appears to chop.\n"); |
---|
830 | } else |
---|
831 | printf ("* is neither chopped nor correctly rounded.\n"); |
---|
832 | if ((RMult == Rounded) && (GMult == No)) |
---|
833 | notify ("Multiplication"); |
---|
834 | } else |
---|
835 | printf ("* is neither chopped nor correctly rounded.\n"); |
---|
836 | /*=============================================*/ |
---|
837 | Milestone = 45; |
---|
838 | /*=============================================*/ |
---|
839 | Y2 = One + U2; |
---|
840 | Y1 = One - U2; |
---|
841 | Z = OneAndHalf + U2 + U2; |
---|
842 | X = Z / Y2; |
---|
843 | T = OneAndHalf - U2 - U2; |
---|
844 | Y = (T - U2) / Y1; |
---|
845 | Z = (Z + U2) / Y2; |
---|
846 | X = X - OneAndHalf; |
---|
847 | Y = Y - T; |
---|
848 | T = T / Y1; |
---|
849 | Z = Z - (OneAndHalf + U2); |
---|
850 | T = (U2 - OneAndHalf) + T; |
---|
851 | if (!((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) { |
---|
852 | X = OneAndHalf / Y2; |
---|
853 | Y = OneAndHalf - U2; |
---|
854 | Z = OneAndHalf + U2; |
---|
855 | X = X - Y; |
---|
856 | T = OneAndHalf / Y1; |
---|
857 | Y = Y / Y1; |
---|
858 | T = T - (Z + U2); |
---|
859 | Y = Y - Z; |
---|
860 | Z = Z / Y2; |
---|
861 | Y1 = (Y2 + U2) / Y2; |
---|
862 | Z = Z - OneAndHalf; |
---|
863 | Y2 = Y1 - Y2; |
---|
864 | Y1 = (F9 - U1) / F9; |
---|
865 | if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero) |
---|
866 | && (Y2 == Zero) && (Y2 == Zero) |
---|
867 | && (Y1 - Half == F9 - Half)) { |
---|
868 | RDiv = Rounded; |
---|
869 | printf ("Division appears to round correctly.\n"); |
---|
870 | if (GDiv == No) |
---|
871 | notify ("Division"); |
---|
872 | } else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero) |
---|
873 | && (Y2 < Zero) && (Y1 - Half < F9 - Half)) { |
---|
874 | RDiv = Chopped; |
---|
875 | printf ("Division appears to chop.\n"); |
---|
876 | } |
---|
877 | } |
---|
878 | if (RDiv == Other) |
---|
879 | printf ("/ is neither chopped nor correctly rounded.\n"); |
---|
880 | BInvrse = One / Radix; |
---|
881 | TstCond (Failure, (BInvrse * Radix - Half == Half), |
---|
882 | "Radix * ( 1 / Radix ) differs from 1"); |
---|
883 | /*=============================================*/ |
---|
884 | Milestone = 50; |
---|
885 | /*=============================================*/ |
---|
886 | TstCond (Failure, ((F9 + U1) - Half == Half) |
---|
887 | && ((BMinusU2 + U2) - One == Radix - One), |
---|
888 | "Incomplete carry-propagation in Addition"); |
---|
889 | X = One - U1 * U1; |
---|
890 | Y = One + U2 * (One - U2); |
---|
891 | Z = F9 - Half; |
---|
892 | X = (X - Half) - Z; |
---|
893 | Y = Y - One; |
---|
894 | if ((X == Zero) && (Y == Zero)) { |
---|
895 | RAddSub = Chopped; |
---|
896 | printf ("Add/Subtract appears to be chopped.\n"); |
---|
897 | } |
---|
898 | if (GAddSub == Yes) { |
---|
899 | X = (Half + U2) * U2; |
---|
900 | Y = (Half - U2) * U2; |
---|
901 | X = One + X; |
---|
902 | Y = One + Y; |
---|
903 | X = (One + U2) - X; |
---|
904 | Y = One - Y; |
---|
905 | if ((X == Zero) && (Y == Zero)) { |
---|
906 | X = (Half + U2) * U1; |
---|
907 | Y = (Half - U2) * U1; |
---|
908 | X = One - X; |
---|
909 | Y = One - Y; |
---|
910 | X = F9 - X; |
---|
911 | Y = One - Y; |
---|
912 | if ((X == Zero) && (Y == Zero)) { |
---|
913 | RAddSub = Rounded; |
---|
914 | printf ("Addition/Subtraction appears to round correctly.\n"); |
---|
915 | if (GAddSub == No) |
---|
916 | notify ("Add/Subtract"); |
---|
917 | } else |
---|
918 | printf ("Addition/Subtraction neither rounds nor chops.\n"); |
---|
919 | } else |
---|
920 | printf ("Addition/Subtraction neither rounds nor chops.\n"); |
---|
921 | } else |
---|
922 | printf ("Addition/Subtraction neither rounds nor chops.\n"); |
---|
923 | S = One; |
---|
924 | X = One + Half * (One + Half); |
---|
925 | Y = (One + U2) * Half; |
---|
926 | Z = X - Y; |
---|
927 | T = Y - X; |
---|
928 | StickyBit = Z + T; |
---|
929 | if (StickyBit != Zero) { |
---|
930 | S = Zero; |
---|
931 | BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n"); |
---|
932 | } |
---|
933 | StickyBit = Zero; |
---|
934 | if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes) |
---|
935 | && (RMult == Rounded) && (RDiv == Rounded) |
---|
936 | && (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2)) { |
---|
937 | printf ("Checking for sticky bit.\n"); |
---|
938 | X = (Half + U1) * U2; |
---|
939 | Y = Half * U2; |
---|
940 | Z = One + Y; |
---|
941 | T = One + X; |
---|
942 | if ((Z - One <= Zero) && (T - One >= U2)) { |
---|
943 | Z = T + Y; |
---|
944 | Y = Z - X; |
---|
945 | if ((Z - T >= U2) && (Y - T == Zero)) { |
---|
946 | X = (Half + U1) * U1; |
---|
947 | Y = Half * U1; |
---|
948 | Z = One - Y; |
---|
949 | T = One - X; |
---|
950 | if ((Z - One == Zero) && (T - F9 == Zero)) { |
---|
951 | Z = (Half - U1) * U1; |
---|
952 | T = F9 - Z; |
---|
953 | Q = F9 - Y; |
---|
954 | if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) { |
---|
955 | Z = (One + U2) * OneAndHalf; |
---|
956 | T = (OneAndHalf + U2) - Z + U2; |
---|
957 | X = One + Half / Radix; |
---|
958 | Y = One + Radix * U2; |
---|
959 | Z = X * Y; |
---|
960 | if (T == Zero && X + Radix * U2 - Z == Zero) { |
---|
961 | if (Radix != Two) { |
---|
962 | X = Two + U2; |
---|
963 | Y = X / Two; |
---|
964 | if ((Y - One == Zero)) |
---|
965 | StickyBit = S; |
---|
966 | } else |
---|
967 | StickyBit = S; |
---|
968 | } |
---|
969 | } |
---|
970 | } |
---|
971 | } |
---|
972 | } |
---|
973 | } |
---|
974 | if (StickyBit == One) |
---|
975 | printf ("Sticky bit apparently used correctly.\n"); |
---|
976 | else |
---|
977 | printf ("Sticky bit used incorrectly or not at all.\n"); |
---|
978 | TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No || |
---|
979 | RMult == Other || RDiv == Other || RAddSub == Other), |
---|
980 | "lack(s) of guard digits or failure(s) to correctly round or chop\n\ |
---|
981 | (noted above) count as one flaw in the final tally below"); |
---|
982 | /*=============================================*/ |
---|
983 | Milestone = 60; |
---|
984 | /*=============================================*/ |
---|
985 | printf ("\n"); |
---|
986 | printf ("Does Multiplication commute? "); |
---|
987 | printf ("Testing on %d random pairs.\n", NoTrials); |
---|
988 | Random9 = SQRT (3.0); |
---|
989 | Random1 = Third; |
---|
990 | I = 1; |
---|
991 | do { |
---|
992 | X = Random (); |
---|
993 | Y = Random (); |
---|
994 | Z9 = Y * X; |
---|
995 | Z = X * Y; |
---|
996 | Z9 = Z - Z9; |
---|
997 | I = I + 1; |
---|
998 | } |
---|
999 | while (!((I > NoTrials) || (Z9 != Zero))); |
---|
1000 | if (I == NoTrials) { |
---|
1001 | Random1 = One + Half / Three; |
---|
1002 | Random2 = (U2 + U1) + One; |
---|
1003 | Z = Random1 * Random2; |
---|
1004 | Y = Random2 * Random1; |
---|
1005 | Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half / |
---|
1006 | Three) * ((U2 + U1) + One); |
---|
1007 | } |
---|
1008 | if (!((I == NoTrials) || (Z9 == Zero))) |
---|
1009 | BadCond (Defect, "X * Y == Y * X trial fails.\n"); |
---|
1010 | else |
---|
1011 | printf (" No failures found in %d integer pairs.\n", NoTrials); |
---|
1012 | /*=============================================*/ |
---|
1013 | Milestone = 70; |
---|
1014 | /*=============================================*/ |
---|
1015 | printf ("\nRunning test of square root(x).\n"); |
---|
1016 | TstCond (Failure, (Zero == SQRT (Zero)) |
---|
1017 | && (-Zero == SQRT (-Zero)) |
---|
1018 | && (One == SQRT (One)), "Square root of 0.0, -0.0 or 1.0 wrong"); |
---|
1019 | MinSqEr = Zero; |
---|
1020 | MaxSqEr = Zero; |
---|
1021 | J = Zero; |
---|
1022 | X = Radix; |
---|
1023 | OneUlp = U2; |
---|
1024 | SqXMinX (Serious); |
---|
1025 | X = BInvrse; |
---|
1026 | OneUlp = BInvrse * U1; |
---|
1027 | SqXMinX (Serious); |
---|
1028 | X = U1; |
---|
1029 | OneUlp = U1 * U1; |
---|
1030 | SqXMinX (Serious); |
---|
1031 | if (J != Zero) |
---|
1032 | Pause (); |
---|
1033 | printf ("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials); |
---|
1034 | J = Zero; |
---|
1035 | X = Two; |
---|
1036 | Y = Radix; |
---|
1037 | if ((Radix != One)) |
---|
1038 | do { |
---|
1039 | X = Y; |
---|
1040 | Y = Radix * Y; |
---|
1041 | } |
---|
1042 | while (!((Y - X >= NoTrials))); |
---|
1043 | OneUlp = X * U2; |
---|
1044 | I = 1; |
---|
1045 | while (I <= NoTrials) { |
---|
1046 | X = X + One; |
---|
1047 | SqXMinX (Defect); |
---|
1048 | if (J > Zero) |
---|
1049 | break; |
---|
1050 | I = I + 1; |
---|
1051 | } |
---|
1052 | printf ("Test for sqrt monotonicity.\n"); |
---|
1053 | I = -1; |
---|
1054 | X = BMinusU2; |
---|
1055 | Y = Radix; |
---|
1056 | Z = Radix + Radix * U2; |
---|
1057 | NotMonot = False; |
---|
1058 | Monot = False; |
---|
1059 | while (!(NotMonot || Monot)) { |
---|
1060 | I = I + 1; |
---|
1061 | X = SQRT (X); |
---|
1062 | Q = SQRT (Y); |
---|
1063 | Z = SQRT (Z); |
---|
1064 | if ((X > Q) || (Q > Z)) |
---|
1065 | NotMonot = True; |
---|
1066 | else { |
---|
1067 | Q = FLOOR (Q + Half); |
---|
1068 | if ((I > 0) || (Radix == Q * Q)) |
---|
1069 | Monot = True; |
---|
1070 | else if (I > 0) { |
---|
1071 | if (I > 1) |
---|
1072 | Monot = True; |
---|
1073 | else { |
---|
1074 | Y = Y * BInvrse; |
---|
1075 | X = Y - U1; |
---|
1076 | Z = Y + U1; |
---|
1077 | } |
---|
1078 | } else { |
---|
1079 | Y = Q; |
---|
1080 | X = Y - U2; |
---|
1081 | Z = Y + U2; |
---|
1082 | } |
---|
1083 | } |
---|
1084 | } |
---|
1085 | if (Monot) |
---|
1086 | printf ("sqrt has passed a test for Monotonicity.\n"); |
---|
1087 | else { |
---|
1088 | BadCond (Defect, ""); |
---|
1089 | printf ("sqrt(X) is non-monotonic for X near %.7e .\n", Y); |
---|
1090 | } |
---|
1091 | /*=============================================*/ |
---|
1092 | Milestone = 80; |
---|
1093 | /*=============================================*/ |
---|
1094 | MinSqEr = MinSqEr + Half; |
---|
1095 | MaxSqEr = MaxSqEr - Half; |
---|
1096 | Y = (SQRT (One + U2) - One) / U2; |
---|
1097 | SqEr = (Y - One) + U2 / Eight; |
---|
1098 | if (SqEr > MaxSqEr) |
---|
1099 | MaxSqEr = SqEr; |
---|
1100 | SqEr = Y + U2 / Eight; |
---|
1101 | if (SqEr < MinSqEr) |
---|
1102 | MinSqEr = SqEr; |
---|
1103 | Y = ((SQRT (F9) - U2) - (One - U2)) / U1; |
---|
1104 | SqEr = Y + U1 / Eight; |
---|
1105 | if (SqEr > MaxSqEr) |
---|
1106 | MaxSqEr = SqEr; |
---|
1107 | SqEr = (Y + One) + U1 / Eight; |
---|
1108 | if (SqEr < MinSqEr) |
---|
1109 | MinSqEr = SqEr; |
---|
1110 | OneUlp = U2; |
---|
1111 | X = OneUlp; |
---|
1112 | for (Indx = 1; Indx <= 3; ++Indx) { |
---|
1113 | Y = SQRT ((X + U1 + X) + F9); |
---|
1114 | Y = ((Y - U2) - ((One - U2) + X)) / OneUlp; |
---|
1115 | Z = ((U1 - X) + F9) * Half * X * X / OneUlp; |
---|
1116 | SqEr = (Y + Half) + Z; |
---|
1117 | if (SqEr < MinSqEr) |
---|
1118 | MinSqEr = SqEr; |
---|
1119 | SqEr = (Y - Half) + Z; |
---|
1120 | if (SqEr > MaxSqEr) |
---|
1121 | MaxSqEr = SqEr; |
---|
1122 | if (((Indx == 1) || (Indx == 3))) |
---|
1123 | X = OneUlp * Sign (X) * FLOOR (Eight / (Nine * SQRT (OneUlp))); |
---|
1124 | else { |
---|
1125 | OneUlp = U1; |
---|
1126 | X = -OneUlp; |
---|
1127 | } |
---|
1128 | } |
---|
1129 | /*=============================================*/ |
---|
1130 | Milestone = 85; |
---|
1131 | /*=============================================*/ |
---|
1132 | SqRWrng = False; |
---|
1133 | Anomaly = False; |
---|
1134 | RSqrt = Other; /* ~dgh */ |
---|
1135 | if (Radix != One) { |
---|
1136 | printf ("Testing whether sqrt is rounded or chopped.\n"); |
---|
1137 | D = FLOOR (Half + POW (Radix, One + Precision - FLOOR (Precision))); |
---|
1138 | /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */ |
---|
1139 | X = D / Radix; |
---|
1140 | Y = D / A1; |
---|
1141 | if ((X != FLOOR (X)) || (Y != FLOOR (Y))) { |
---|
1142 | Anomaly = True; |
---|
1143 | } else { |
---|
1144 | X = Zero; |
---|
1145 | Z2 = X; |
---|
1146 | Y = One; |
---|
1147 | Y2 = Y; |
---|
1148 | Z1 = Radix - One; |
---|
1149 | FourD = Four * D; |
---|
1150 | do { |
---|
1151 | if (Y2 > Z2) { |
---|
1152 | Q = Radix; |
---|
1153 | Y1 = Y; |
---|
1154 | do { |
---|
1155 | X1 = FABS (Q + FLOOR (Half - Q / Y1) * Y1); |
---|
1156 | Q = Y1; |
---|
1157 | Y1 = X1; |
---|
1158 | } |
---|
1159 | while (!(X1 <= Zero)); |
---|
1160 | if (Q <= One) { |
---|
1161 | Z2 = Y2; |
---|
1162 | Z = Y; |
---|
1163 | } |
---|
1164 | } |
---|
1165 | Y = Y + Two; |
---|
1166 | X = X + Eight; |
---|
1167 | Y2 = Y2 + X; |
---|
1168 | if (Y2 >= FourD) |
---|
1169 | Y2 = Y2 - FourD; |
---|
1170 | } |
---|
1171 | while (!(Y >= D)); |
---|
1172 | X8 = FourD - Z2; |
---|
1173 | Q = (X8 + Z * Z) / FourD; |
---|
1174 | X8 = X8 / Eight; |
---|
1175 | if (Q != FLOOR (Q)) |
---|
1176 | Anomaly = True; |
---|
1177 | else { |
---|
1178 | Break = False; |
---|
1179 | do { |
---|
1180 | X = Z1 * Z; |
---|
1181 | X = X - FLOOR (X / Radix) * Radix; |
---|
1182 | if (X == One) |
---|
1183 | Break = True; |
---|
1184 | else |
---|
1185 | Z1 = Z1 - One; |
---|
1186 | } |
---|
1187 | while (!(Break || (Z1 <= Zero))); |
---|
1188 | if ((Z1 <= Zero) && (!Break)) |
---|
1189 | Anomaly = True; |
---|
1190 | else { |
---|
1191 | if (Z1 > RadixD2) |
---|
1192 | Z1 = Z1 - Radix; |
---|
1193 | do { |
---|
1194 | NewD (); |
---|
1195 | } |
---|
1196 | while (!(U2 * D >= F9)); |
---|
1197 | if (D * Radix - D != WVar - D) |
---|
1198 | Anomaly = True; |
---|
1199 | else { |
---|
1200 | Z2 = D; |
---|
1201 | I = 0; |
---|
1202 | Y = D + (One + Z) * Half; |
---|
1203 | X = D + Z + Q; |
---|
1204 | SR3750 (); |
---|
1205 | Y = D + (One - Z) * Half + D; |
---|
1206 | X = D - Z + D; |
---|
1207 | X = X + Q + X; |
---|
1208 | SR3750 (); |
---|
1209 | NewD (); |
---|
1210 | if (D - Z2 != WVar - Z2) |
---|
1211 | Anomaly = True; |
---|
1212 | else { |
---|
1213 | Y = (D - Z2) + (Z2 + (One - Z) * Half); |
---|
1214 | X = (D - Z2) + (Z2 - Z + Q); |
---|
1215 | SR3750 (); |
---|
1216 | Y = (One + Z) * Half; |
---|
1217 | X = Q; |
---|
1218 | SR3750 (); |
---|
1219 | if (I == 0) |
---|
1220 | Anomaly = True; |
---|
1221 | } |
---|
1222 | } |
---|
1223 | } |
---|
1224 | } |
---|
1225 | } |
---|
1226 | if ((I == 0) || Anomaly) { |
---|
1227 | BadCond (Failure, "Anomalous arithmetic with Integer < "); |
---|
1228 | printf ("Radix^Precision = %.7e\n", WVar); |
---|
1229 | printf (" fails test whether sqrt rounds or chops.\n"); |
---|
1230 | SqRWrng = True; |
---|
1231 | } |
---|
1232 | } |
---|
1233 | if (!Anomaly) { |
---|
1234 | if (!((MinSqEr < Zero) || (MaxSqEr > Zero))) { |
---|
1235 | RSqrt = Rounded; |
---|
1236 | printf ("Square root appears to be correctly rounded.\n"); |
---|
1237 | } else { |
---|
1238 | if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half) |
---|
1239 | || (MinSqEr + Radix < Half)) |
---|
1240 | SqRWrng = True; |
---|
1241 | else { |
---|
1242 | RSqrt = Chopped; |
---|
1243 | printf ("Square root appears to be chopped.\n"); |
---|
1244 | } |
---|
1245 | } |
---|
1246 | } |
---|
1247 | if (SqRWrng) { |
---|
1248 | printf ("Square root is neither chopped nor correctly rounded.\n"); |
---|
1249 | printf ("Observed errors run from %.7e ", MinSqEr - Half); |
---|
1250 | printf ("to %.7e ulps.\n", Half + MaxSqEr); |
---|
1251 | TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix, |
---|
1252 | "sqrt gets too many last digits wrong"); |
---|
1253 | } |
---|
1254 | /*=============================================*/ |
---|
1255 | Milestone = 90; |
---|
1256 | /*=============================================*/ |
---|
1257 | Pause (); |
---|
1258 | printf ("Testing powers Z^i for small Integers Z and i.\n"); |
---|
1259 | N = 0; |
---|
1260 | /* ... test powers of zero. */ |
---|
1261 | I = 0; |
---|
1262 | Z = -Zero; |
---|
1263 | M = 3; |
---|
1264 | Break = False; |
---|
1265 | do { |
---|
1266 | X = One; |
---|
1267 | SR3980 (); |
---|
1268 | if (I <= 10) { |
---|
1269 | I = 1023; |
---|
1270 | SR3980 (); |
---|
1271 | } |
---|
1272 | if (Z == MinusOne) |
---|
1273 | Break = True; |
---|
1274 | else { |
---|
1275 | Z = MinusOne; |
---|
1276 | /* .. if(-1)^N is invalid, replace MinusOne by One. */ |
---|
1277 | I = -4; |
---|
1278 | } |
---|
1279 | } |
---|
1280 | while (!Break); |
---|
1281 | PrintIfNPositive (); |
---|
1282 | N1 = N; |
---|
1283 | N = 0; |
---|
1284 | Z = A1; |
---|
1285 | M = (int) FLOOR (Two * LOG (WVar) / LOG (A1)); |
---|
1286 | Break = False; |
---|
1287 | do { |
---|
1288 | X = Z; |
---|
1289 | I = 1; |
---|
1290 | SR3980 (); |
---|
1291 | if (Z == AInvrse) |
---|
1292 | Break = True; |
---|
1293 | else |
---|
1294 | Z = AInvrse; |
---|
1295 | } |
---|
1296 | while (!(Break)); |
---|
1297 | /*=============================================*/ |
---|
1298 | Milestone = 100; |
---|
1299 | /*=============================================*/ |
---|
1300 | /* Powers of Radix have been tested, */ |
---|
1301 | /* next try a few primes */ |
---|
1302 | M = NoTrials; |
---|
1303 | Z = Three; |
---|
1304 | do { |
---|
1305 | X = Z; |
---|
1306 | I = 1; |
---|
1307 | SR3980 (); |
---|
1308 | do { |
---|
1309 | Z = Z + Two; |
---|
1310 | } |
---|
1311 | while (Three * FLOOR (Z / Three) == Z); |
---|
1312 | } |
---|
1313 | while (Z < Eight * Three); |
---|
1314 | if (N > 0) { |
---|
1315 | printf ("Errors like this may invalidate financial calculations\n"); |
---|
1316 | printf ("\tinvolving interest rates.\n"); |
---|
1317 | } |
---|
1318 | PrintIfNPositive (); |
---|
1319 | N += N1; |
---|
1320 | if (N == 0) |
---|
1321 | printf ("... no discrepancies found.\n"); |
---|
1322 | if (N > 0) |
---|
1323 | Pause (); |
---|
1324 | else |
---|
1325 | printf ("\n"); |
---|
1326 | /*=============================================*/ |
---|
1327 | Milestone = 110; |
---|
1328 | /*=============================================*/ |
---|
1329 | printf ("Seeking Underflow thresholds UfThold and E0.\n"); |
---|
1330 | D = U1; |
---|
1331 | if (Precision != FLOOR (Precision)) { |
---|
1332 | D = BInvrse; |
---|
1333 | X = Precision; |
---|
1334 | do { |
---|
1335 | D = D * BInvrse; |
---|
1336 | X = X - One; |
---|
1337 | } |
---|
1338 | while (X > Zero); |
---|
1339 | } |
---|
1340 | Y = One; |
---|
1341 | Z = D; |
---|
1342 | /* ... D is power of 1/Radix < 1. */ |
---|
1343 | do { |
---|
1344 | C = Y; |
---|
1345 | Y = Z; |
---|
1346 | Z = Y * Y; |
---|
1347 | } |
---|
1348 | while ((Y > Z) && (Z + Z > Z)); |
---|
1349 | Y = C; |
---|
1350 | Z = Y * D; |
---|
1351 | do { |
---|
1352 | C = Y; |
---|
1353 | Y = Z; |
---|
1354 | Z = Y * D; |
---|
1355 | } |
---|
1356 | while ((Y > Z) && (Z + Z > Z)); |
---|
1357 | if (Radix < Two) |
---|
1358 | HInvrse = Two; |
---|
1359 | else |
---|
1360 | HInvrse = Radix; |
---|
1361 | HVar = One / HInvrse; |
---|
1362 | /* ... 1/HInvrse == HVar == Min(1/Radix, 1/2) */ |
---|
1363 | CInvrse = One / C; |
---|
1364 | E0 = C; |
---|
1365 | Z = E0 * HVar; |
---|
1366 | /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */ |
---|
1367 | do { |
---|
1368 | Y = E0; |
---|
1369 | E0 = Z; |
---|
1370 | Z = E0 * HVar; |
---|
1371 | } |
---|
1372 | while ((E0 > Z) && (Z + Z > Z)); |
---|
1373 | UfThold = E0; |
---|
1374 | E1 = Zero; |
---|
1375 | Q = Zero; |
---|
1376 | E9 = U2; |
---|
1377 | S = One + E9; |
---|
1378 | D = C * S; |
---|
1379 | if (D <= C) { |
---|
1380 | E9 = Radix * U2; |
---|
1381 | S = One + E9; |
---|
1382 | D = C * S; |
---|
1383 | if (D <= C) { |
---|
1384 | BadCond (Failure, "multiplication gets too many last digits wrong.\n"); |
---|
1385 | Underflow = E0; |
---|
1386 | Y1 = Zero; |
---|
1387 | PseudoZero = Z; |
---|
1388 | Pause (); |
---|
1389 | } |
---|
1390 | } else { |
---|
1391 | Underflow = D; |
---|
1392 | PseudoZero = Underflow * HVar; |
---|
1393 | UfThold = Zero; |
---|
1394 | do { |
---|
1395 | Y1 = Underflow; |
---|
1396 | Underflow = PseudoZero; |
---|
1397 | if (E1 + E1 <= E1) { |
---|
1398 | Y2 = Underflow * HInvrse; |
---|
1399 | E1 = FABS (Y1 - Y2); |
---|
1400 | Q = Y1; |
---|
1401 | if ((UfThold == Zero) && (Y1 != Y2)) |
---|
1402 | UfThold = Y1; |
---|
1403 | } |
---|
1404 | PseudoZero = PseudoZero * HVar; |
---|
1405 | } |
---|
1406 | while ((Underflow > PseudoZero) |
---|
1407 | && (PseudoZero + PseudoZero > PseudoZero)); |
---|
1408 | } |
---|
1409 | /* Comment line 4530 .. 4560 */ |
---|
1410 | if (PseudoZero != Zero) { |
---|
1411 | printf ("\n"); |
---|
1412 | Z = PseudoZero; |
---|
1413 | /* ... Test PseudoZero for "phoney- zero" violates */ |
---|
1414 | /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero |
---|
1415 | ... */ |
---|
1416 | if (PseudoZero <= Zero) { |
---|
1417 | BadCond (Failure, "Positive expressions can underflow to an\n"); |
---|
1418 | printf ("allegedly negative value\n"); |
---|
1419 | printf ("PseudoZero that prints out as: %g .\n", PseudoZero); |
---|
1420 | X = -PseudoZero; |
---|
1421 | if (X <= Zero) { |
---|
1422 | printf ("But -PseudoZero, which should be\n"); |
---|
1423 | printf ("positive, isn't; it prints out as %g .\n", X); |
---|
1424 | } |
---|
1425 | } else { |
---|
1426 | BadCond (Flaw, "Underflow can stick at an allegedly positive\n"); |
---|
1427 | printf ("value PseudoZero that prints out as %g .\n", PseudoZero); |
---|
1428 | } |
---|
1429 | TstPtUf (); |
---|
1430 | } |
---|
1431 | /*=============================================*/ |
---|
1432 | Milestone = 120; |
---|
1433 | /*=============================================*/ |
---|
1434 | if (CInvrse * Y > CInvrse * Y1) { |
---|
1435 | S = HVar * S; |
---|
1436 | E0 = Underflow; |
---|
1437 | } |
---|
1438 | if (!((E1 == Zero) || (E1 == E0))) { |
---|
1439 | BadCond (Defect, ""); |
---|
1440 | if (E1 < E0) { |
---|
1441 | printf ("Products underflow at a higher"); |
---|
1442 | printf (" threshold than differences.\n"); |
---|
1443 | if (PseudoZero == Zero) |
---|
1444 | E0 = E1; |
---|
1445 | } else { |
---|
1446 | printf ("Difference underflows at a higher"); |
---|
1447 | printf (" threshold than products.\n"); |
---|
1448 | } |
---|
1449 | } |
---|
1450 | printf ("Smallest strictly positive number found is E0 = %g .\n", E0); |
---|
1451 | Z = E0; |
---|
1452 | TstPtUf (); |
---|
1453 | Underflow = E0; |
---|
1454 | if (N == 1) |
---|
1455 | Underflow = Y; |
---|
1456 | I = 4; |
---|
1457 | if (E1 == Zero) |
---|
1458 | I = 3; |
---|
1459 | if (UfThold == Zero) |
---|
1460 | I = I - 2; |
---|
1461 | UfNGrad = True; |
---|
1462 | switch (I) { |
---|
1463 | case 1: |
---|
1464 | UfThold = Underflow; |
---|
1465 | if ((CInvrse * Q) != ((CInvrse * Y) * S)) { |
---|
1466 | UfThold = Y; |
---|
1467 | BadCond (Failure, "Either accuracy deteriorates as numbers\n"); |
---|
1468 | printf ("approach a threshold = %.17e\n", UfThold);; |
---|
1469 | printf (" coming down from %.17e\n", C); |
---|
1470 | printf (" or else multiplication gets too many last digits wrong.\n"); |
---|
1471 | } |
---|
1472 | Pause (); |
---|
1473 | break; |
---|
1474 | |
---|
1475 | case 2: |
---|
1476 | BadCond (Failure, "Underflow confuses Comparison, which alleges that\n"); |
---|
1477 | printf ("Q == Y while denying that |Q - Y| == 0; these values\n"); |
---|
1478 | printf ("print out as Q = %.17e, Y = %.17e .\n", Q, Y2); |
---|
1479 | printf ("|Q - Y| = %.17e .\n", FABS (Q - Y2)); |
---|
1480 | UfThold = Q; |
---|
1481 | break; |
---|
1482 | |
---|
1483 | case 3: |
---|
1484 | X = X; |
---|
1485 | break; |
---|
1486 | |
---|
1487 | case 4: |
---|
1488 | if ((Q == UfThold) && (E1 == E0) |
---|
1489 | && (FABS (UfThold - E1 / E9) <= E1)) { |
---|
1490 | UfNGrad = False; |
---|
1491 | printf ("Underflow is gradual; it incurs Absolute Error =\n"); |
---|
1492 | printf ("(roundoff in UfThold) < E0.\n"); |
---|
1493 | Y = E0 * CInvrse; |
---|
1494 | Y = Y * (OneAndHalf + U2); |
---|
1495 | X = CInvrse * (One + U2); |
---|
1496 | Y = Y / X; |
---|
1497 | IEEE = (Y == E0); |
---|
1498 | } |
---|
1499 | } |
---|
1500 | if (UfNGrad) { |
---|
1501 | printf ("\n"); |
---|
1502 | sigsave = _sigfpe; |
---|
1503 | if (setjmp (ovfl_buf)) { |
---|
1504 | printf ("Underflow / UfThold failed!\n"); |
---|
1505 | R = HVar + HVar; |
---|
1506 | } else |
---|
1507 | R = SQRT (Underflow / UfThold); |
---|
1508 | sigsave = 0; |
---|
1509 | if (R <= HVar) { |
---|
1510 | Z = R * UfThold; |
---|
1511 | X = Z * (One + R * HVar * (One + HVar)); |
---|
1512 | } else { |
---|
1513 | Z = UfThold; |
---|
1514 | X = Z * (One + HVar * HVar * (One + HVar)); |
---|
1515 | } |
---|
1516 | if (!((X == Z) || (X - Z != Zero))) { |
---|
1517 | BadCond (Flaw, ""); |
---|
1518 | printf ("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z); |
---|
1519 | Z9 = X - Z; |
---|
1520 | printf ("yet X - Z yields %.17e .\n", Z9); |
---|
1521 | printf (" Should this NOT signal Underflow, "); |
---|
1522 | printf ("this is a SERIOUS DEFECT\nthat causes "); |
---|
1523 | printf ("confusion when innocent statements like\n");; |
---|
1524 | printf (" if (X == Z) ... else"); |
---|
1525 | printf (" ... (f(X) - f(Z)) / (X - Z) ...\n"); |
---|
1526 | printf ("encounter Division by Zero although actually\n"); |
---|
1527 | sigsave = _sigfpe; |
---|
1528 | if (setjmp (ovfl_buf)) |
---|
1529 | printf ("X / Z fails!\n"); |
---|
1530 | else |
---|
1531 | printf ("X / Z = 1 + %g .\n", (X / Z - Half) - Half); |
---|
1532 | sigsave = 0; |
---|
1533 | } |
---|
1534 | } |
---|
1535 | printf ("The Underflow threshold is %.17e, %s\n", UfThold, |
---|
1536 | " below which"); |
---|
1537 | printf ("calculation may suffer larger Relative error than "); |
---|
1538 | printf ("merely roundoff.\n"); |
---|
1539 | Y2 = U1 * U1; |
---|
1540 | Y = Y2 * Y2; |
---|
1541 | Y2 = Y * U1; |
---|
1542 | if (Y2 <= UfThold) { |
---|
1543 | if (Y > E0) { |
---|
1544 | BadCond (Defect, ""); |
---|
1545 | I = 5; |
---|
1546 | } else { |
---|
1547 | BadCond (Serious, ""); |
---|
1548 | I = 4; |
---|
1549 | } |
---|
1550 | printf ("Range is too narrow; U1^%d Underflows.\n", I); |
---|
1551 | } |
---|
1552 | /*=============================================*/ |
---|
1553 | Milestone = 130; |
---|
1554 | /*=============================================*/ |
---|
1555 | Y = -FLOOR (Half - TwoForty * LOG (UfThold) / LOG (HInvrse)) / TwoForty; |
---|
1556 | Y2 = Y + Y; |
---|
1557 | printf ("Since underflow occurs below the threshold\n"); |
---|
1558 | printf ("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y); |
---|
1559 | printf ("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", |
---|
1560 | HInvrse, Y2); |
---|
1561 | printf ("actually calculating yields:"); |
---|
1562 | if (setjmp (ovfl_buf)) { |
---|
1563 | sigsave = 0; |
---|
1564 | BadCond (Serious, "trap on underflow.\n"); |
---|
1565 | } else { |
---|
1566 | sigsave = _sigfpe; |
---|
1567 | V9 = POW (HInvrse, Y2); |
---|
1568 | sigsave = 0; |
---|
1569 | printf (" %.17e .\n", V9); |
---|
1570 | if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) { |
---|
1571 | BadCond (Serious, "this is not between 0 and underflow\n"); |
---|
1572 | printf (" threshold = %.17e .\n", UfThold); |
---|
1573 | } else if (!(V9 > UfThold * (One + E9))) |
---|
1574 | printf ("This computed value is O.K.\n"); |
---|
1575 | else { |
---|
1576 | BadCond (Defect, "this is not between 0 and underflow\n"); |
---|
1577 | printf (" threshold = %.17e .\n", UfThold); |
---|
1578 | } |
---|
1579 | } |
---|
1580 | /*=============================================*/ |
---|
1581 | Milestone = 140; |
---|
1582 | /*=============================================*/ |
---|
1583 | printf ("\n"); |
---|
1584 | /* ...calculate Exp2 == exp(2) == 7.389056099... */ |
---|
1585 | X = Zero; |
---|
1586 | I = 2; |
---|
1587 | Y = Two * Three; |
---|
1588 | Q = Zero; |
---|
1589 | N = 0; |
---|
1590 | do { |
---|
1591 | Z = X; |
---|
1592 | I = I + 1; |
---|
1593 | Y = Y / (I + I); |
---|
1594 | R = Y + Q; |
---|
1595 | X = Z + R; |
---|
1596 | Q = (Z - X) + R; |
---|
1597 | } |
---|
1598 | while (X > Z); |
---|
1599 | Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo); |
---|
1600 | X = Z * Z; |
---|
1601 | Exp2 = X * X; |
---|
1602 | X = F9; |
---|
1603 | Y = X - U1; |
---|
1604 | printf ("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n", |
---|
1605 | Exp2); |
---|
1606 | for (I = 1;;) { |
---|
1607 | Z = X - BInvrse; |
---|
1608 | Z = (X + One) / (Z - (One - BInvrse)); |
---|
1609 | Q = POW (X, Z) - Exp2; |
---|
1610 | if (FABS (Q) > TwoForty * U2) { |
---|
1611 | N = 1; |
---|
1612 | V9 = (X - BInvrse) - (One - BInvrse); |
---|
1613 | BadCond (Defect, "Calculated"); |
---|
1614 | printf (" %.17e for\n", POW (X, Z)); |
---|
1615 | printf ("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z); |
---|
1616 | printf ("\tdiffers from correct value by %.17e .\n", Q); |
---|
1617 | printf ("\tThis much error may spoil financial\n"); |
---|
1618 | printf ("\tcalculations involving tiny interest rates.\n"); |
---|
1619 | break; |
---|
1620 | } else { |
---|
1621 | Z = (Y - X) * Two + Y; |
---|
1622 | X = Y; |
---|
1623 | Y = Z; |
---|
1624 | Z = One + (X - F9) * (X - F9); |
---|
1625 | if (Z > One && I < NoTrials) |
---|
1626 | I++; |
---|
1627 | else { |
---|
1628 | if (X > One) { |
---|
1629 | if (N == 0) |
---|
1630 | printf ("Accuracy seems adequate.\n"); |
---|
1631 | break; |
---|
1632 | } else { |
---|
1633 | X = One + U2; |
---|
1634 | Y = U2 + U2; |
---|
1635 | Y += X; |
---|
1636 | I = 1; |
---|
1637 | } |
---|
1638 | } |
---|
1639 | } |
---|
1640 | } |
---|
1641 | /*=============================================*/ |
---|
1642 | Milestone = 150; |
---|
1643 | /*=============================================*/ |
---|
1644 | printf ("Testing powers Z^Q at four nearly extreme values.\n"); |
---|
1645 | N = 0; |
---|
1646 | Z = A1; |
---|
1647 | Q = FLOOR (Half - LOG (C) / LOG (A1)); |
---|
1648 | Break = False; |
---|
1649 | do { |
---|
1650 | X = CInvrse; |
---|
1651 | Y = POW (Z, Q); |
---|
1652 | IsYeqX (); |
---|
1653 | Q = -Q; |
---|
1654 | X = C; |
---|
1655 | Y = POW (Z, Q); |
---|
1656 | IsYeqX (); |
---|
1657 | if (Z < One) |
---|
1658 | Break = True; |
---|
1659 | else |
---|
1660 | Z = AInvrse; |
---|
1661 | } |
---|
1662 | while (!(Break)); |
---|
1663 | PrintIfNPositive (); |
---|
1664 | if (N == 0) |
---|
1665 | printf (" ... no discrepancies found.\n"); |
---|
1666 | printf ("\n"); |
---|
1667 | |
---|
1668 | /*=============================================*/ |
---|
1669 | Milestone = 160; |
---|
1670 | /*=============================================*/ |
---|
1671 | Pause (); |
---|
1672 | printf ("Searching for Overflow threshold:\n"); |
---|
1673 | printf ("This may generate an error.\n"); |
---|
1674 | Y = -CInvrse; |
---|
1675 | V9 = HInvrse * Y; |
---|
1676 | sigsave = _sigfpe; |
---|
1677 | if (setjmp (ovfl_buf)) { |
---|
1678 | I = 0; |
---|
1679 | V9 = Y; |
---|
1680 | goto overflow; |
---|
1681 | } |
---|
1682 | do { |
---|
1683 | V = Y; |
---|
1684 | Y = V9; |
---|
1685 | V9 = HInvrse * Y; |
---|
1686 | } |
---|
1687 | while (V9 < Y); |
---|
1688 | I = 1; |
---|
1689 | overflow: |
---|
1690 | sigsave = 0; |
---|
1691 | Z = V9; |
---|
1692 | printf ("Can `Z = -Y' overflow?\n"); |
---|
1693 | printf ("Trying it on Y = %.17e .\n", Y); |
---|
1694 | V9 = -Y; |
---|
1695 | V0 = V9; |
---|
1696 | if (V - Y == V + V0) |
---|
1697 | printf ("Seems O.K.\n"); |
---|
1698 | else { |
---|
1699 | printf ("finds a "); |
---|
1700 | BadCond (Flaw, "-(-Y) differs from Y.\n"); |
---|
1701 | } |
---|
1702 | if (Z != Y) { |
---|
1703 | BadCond (Serious, ""); |
---|
1704 | printf ("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z); |
---|
1705 | } |
---|
1706 | if (I) { |
---|
1707 | Y = V * (HInvrse * U2 - HInvrse); |
---|
1708 | Z = Y + ((One - HInvrse) * U2) * V; |
---|
1709 | if (Z < V0) |
---|
1710 | Y = Z; |
---|
1711 | if (Y < V0) |
---|
1712 | V = Y; |
---|
1713 | if (V0 - V < V0) |
---|
1714 | V = V0; |
---|
1715 | } else { |
---|
1716 | V = Y * (HInvrse * U2 - HInvrse); |
---|
1717 | V = V + ((One - HInvrse) * U2) * Y; |
---|
1718 | } |
---|
1719 | printf ("Overflow threshold is V = %.17e .\n", V); |
---|
1720 | if (I) |
---|
1721 | printf ("Overflow saturates at V0 = %.17e .\n", V0); |
---|
1722 | else |
---|
1723 | printf ("There is no saturation value because \ |
---|
1724 | the system traps on overflow.\n"); |
---|
1725 | V9 = V * One; |
---|
1726 | printf ("No Overflow should be signaled for V * 1 = %.17e\n", V9); |
---|
1727 | V9 = V / One; |
---|
1728 | printf (" nor for V / 1 = %.17e .\n", V9); |
---|
1729 | printf ("Any overflow signal separating this * from the one\n"); |
---|
1730 | printf ("above is a DEFECT.\n"); |
---|
1731 | /*=============================================*/ |
---|
1732 | Milestone = 170; |
---|
1733 | /*=============================================*/ |
---|
1734 | if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) { |
---|
1735 | BadCond (Failure, "Comparisons involving "); |
---|
1736 | printf ("+-%g, +-%g\nand +-%g are confused by Overflow.", |
---|
1737 | V, V0, UfThold); |
---|
1738 | } |
---|
1739 | /*=============================================*/ |
---|
1740 | Milestone = 175; |
---|
1741 | /*=============================================*/ |
---|
1742 | printf ("\n"); |
---|
1743 | for (Indx = 1; Indx <= 3; ++Indx) { |
---|
1744 | switch (Indx) { |
---|
1745 | case 1: |
---|
1746 | Z = UfThold; |
---|
1747 | break; |
---|
1748 | case 2: |
---|
1749 | Z = E0; |
---|
1750 | break; |
---|
1751 | case 3: |
---|
1752 | Z = PseudoZero; |
---|
1753 | break; |
---|
1754 | } |
---|
1755 | if (Z != Zero) { |
---|
1756 | V9 = SQRT (Z); |
---|
1757 | Y = V9 * V9; |
---|
1758 | if (Y / (One - Radix * E9) < Z |
---|
1759 | || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */ |
---|
1760 | if (V9 > U1) |
---|
1761 | BadCond (Serious, ""); |
---|
1762 | else |
---|
1763 | BadCond (Defect, ""); |
---|
1764 | printf ("Comparison alleges that what prints as Z = %.17e\n", Z); |
---|
1765 | printf (" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y); |
---|
1766 | } |
---|
1767 | } |
---|
1768 | } |
---|
1769 | /*=============================================*/ |
---|
1770 | Milestone = 180; |
---|
1771 | /*=============================================*/ |
---|
1772 | for (Indx = 1; Indx <= 2; ++Indx) { |
---|
1773 | if (Indx == 1) |
---|
1774 | Z = V; |
---|
1775 | else |
---|
1776 | Z = V0; |
---|
1777 | V9 = SQRT (Z); |
---|
1778 | X = (One - Radix * E9) * V9; |
---|
1779 | V9 = V9 * X; |
---|
1780 | if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) { |
---|
1781 | Y = V9; |
---|
1782 | if (X < WVar) |
---|
1783 | BadCond (Serious, ""); |
---|
1784 | else |
---|
1785 | BadCond (Defect, ""); |
---|
1786 | printf ("Comparison alleges that Z = %17e\n", Z); |
---|
1787 | printf (" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y); |
---|
1788 | } |
---|
1789 | } |
---|
1790 | /*=============================================*/ |
---|
1791 | Milestone = 190; |
---|
1792 | /*=============================================*/ |
---|
1793 | Pause (); |
---|
1794 | X = UfThold * V; |
---|
1795 | Y = Radix * Radix; |
---|
1796 | if (X * Y < One || X > Y) { |
---|
1797 | if (X * Y < U1 || X > Y / U1) |
---|
1798 | BadCond (Defect, "Badly"); |
---|
1799 | else |
---|
1800 | BadCond (Flaw, ""); |
---|
1801 | |
---|
1802 | printf (" unbalanced range; UfThold * V = %.17e\n\t%s\n", |
---|
1803 | X, "is too far from 1.\n"); |
---|
1804 | } |
---|
1805 | /*=============================================*/ |
---|
1806 | Milestone = 200; |
---|
1807 | /*=============================================*/ |
---|
1808 | for (Indx = 1; Indx <= 5; ++Indx) { |
---|
1809 | X = F9; |
---|
1810 | switch (Indx) { |
---|
1811 | case 2: |
---|
1812 | X = One + U2; |
---|
1813 | break; |
---|
1814 | case 3: |
---|
1815 | X = V; |
---|
1816 | break; |
---|
1817 | case 4: |
---|
1818 | X = UfThold; |
---|
1819 | break; |
---|
1820 | case 5: |
---|
1821 | X = Radix; |
---|
1822 | } |
---|
1823 | Y = X; |
---|
1824 | sigsave = _sigfpe; |
---|
1825 | if (setjmp (ovfl_buf)) |
---|
1826 | printf (" X / X traps when X = %g\n", X); |
---|
1827 | else { |
---|
1828 | V9 = (Y / X - Half) - Half; |
---|
1829 | if (V9 == Zero) |
---|
1830 | continue; |
---|
1831 | if (V9 == -U1 && Indx < 5) |
---|
1832 | BadCond (Flaw, ""); |
---|
1833 | else |
---|
1834 | BadCond (Serious, ""); |
---|
1835 | printf (" X / X differs from 1 when X = %.17e\n", X); |
---|
1836 | printf (" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9); |
---|
1837 | } |
---|
1838 | sigsave = 0; |
---|
1839 | } |
---|
1840 | /*=============================================*/ |
---|
1841 | Milestone = 210; |
---|
1842 | /*=============================================*/ |
---|
1843 | MyZero = Zero; |
---|
1844 | printf ("\n"); |
---|
1845 | printf ("What message and/or values does Division by Zero produce?\n"); |
---|
1846 | #ifndef BATCHMODE |
---|
1847 | printf ("This can interupt your program. You can "); |
---|
1848 | printf ("skip this part if you wish.\n"); |
---|
1849 | printf ("Do you wish to compute 1 / 0? "); |
---|
1850 | fflush (stdout); |
---|
1851 | read (KEYBOARD, ch, 8); |
---|
1852 | if ((ch[0] == 'Y') || (ch[0] == 'y')) { |
---|
1853 | #endif /* !BATCHMODE */ |
---|
1854 | sigsave = _sigfpe; |
---|
1855 | printf (" Trying to compute 1 / 0 produces ..."); |
---|
1856 | if (!setjmp (ovfl_buf)) |
---|
1857 | printf (" %.7e .\n", One / MyZero); |
---|
1858 | sigsave = 0; |
---|
1859 | #ifndef BATCHMODE |
---|
1860 | } else |
---|
1861 | printf ("O.K.\n"); |
---|
1862 | printf ("\nDo you wish to compute 0 / 0? "); |
---|
1863 | fflush (stdout); |
---|
1864 | read (KEYBOARD, ch, 80); |
---|
1865 | if ((ch[0] == 'Y') || (ch[0] == 'y')) { |
---|
1866 | #endif /* !BATCHMODE */ |
---|
1867 | sigsave = _sigfpe; |
---|
1868 | printf ("\n Trying to compute 0 / 0 produces ..."); |
---|
1869 | if (!setjmp (ovfl_buf)) |
---|
1870 | printf (" %.7e .\n", Zero / MyZero); |
---|
1871 | sigsave = 0; |
---|
1872 | #ifndef BATCHMODE |
---|
1873 | } else |
---|
1874 | printf ("O.K.\n"); |
---|
1875 | #endif /* !BATCHMODE */ |
---|
1876 | /*=============================================*/ |
---|
1877 | Milestone = 220; |
---|
1878 | /*=============================================*/ |
---|
1879 | |
---|
1880 | Pause (); |
---|
1881 | printf ("\n"); |
---|
1882 | { |
---|
1883 | static char *msg[] = |
---|
1884 | { |
---|
1885 | "FAILUREs encountered =", |
---|
1886 | "SERIOUS DEFECTs discovered =", |
---|
1887 | "DEFECTs discovered =", |
---|
1888 | "FLAWs discovered ="}; |
---|
1889 | int i; |
---|
1890 | for (i = 0; i < 4; i++) |
---|
1891 | if (ErrCnt[i]) |
---|
1892 | printf ("The number of %-29s %d.\n", |
---|
1893 | msg[i], ErrCnt[i]); |
---|
1894 | } |
---|
1895 | |
---|
1896 | printf ("\n"); |
---|
1897 | if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] |
---|
1898 | + ErrCnt[Flaw]) > 0) { |
---|
1899 | if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[ |
---|
1900 | Defect] == 0) && (ErrCnt[Flaw] > 0)) { |
---|
1901 | printf ("The arithmetic diagnosed seems "); |
---|
1902 | printf ("Satisfactory though flawed.\n"); |
---|
1903 | } |
---|
1904 | if ((ErrCnt[Failure] + ErrCnt[Serious] == 0) |
---|
1905 | && (ErrCnt[Defect] > 0)) { |
---|
1906 | printf ("The arithmetic diagnosed may be Acceptable\n"); |
---|
1907 | printf ("despite inconvenient Defects.\n"); |
---|
1908 | } |
---|
1909 | if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) { |
---|
1910 | printf ("The arithmetic diagnosed has "); |
---|
1911 | printf ("unacceptable Serious Defects.\n"); |
---|
1912 | } |
---|
1913 | if (ErrCnt[Failure] > 0) { |
---|
1914 | printf ("Potentially fatal FAILURE may have spoiled this"); |
---|
1915 | printf (" program's subsequent diagnoses.\n"); |
---|
1916 | } |
---|
1917 | } else { |
---|
1918 | |
---|
1919 | printf ("No failures, defects nor flaws have been discovered.\n"); |
---|
1920 | if (!((RMult == Rounded) && (RDiv == Rounded) |
---|
1921 | && (RAddSub == Rounded) && (RSqrt == Rounded))) |
---|
1922 | printf ("The arithmetic diagnosed seems Satisfactory.\n"); |
---|
1923 | else { |
---|
1924 | if (StickyBit >= One && |
---|
1925 | (Radix - Two) * (Radix - Nine - One) == Zero) { |
---|
1926 | printf ("Rounding appears to conform to "); |
---|
1927 | printf ("the proposed IEEE standard P"); |
---|
1928 | if ((Radix == Two) && |
---|
1929 | ((Precision - Four * Three * Two) * |
---|
1930 | (Precision - TwentySeven - |
---|
1931 | TwentySeven + One) == Zero)) |
---|
1932 | printf ("754"); |
---|
1933 | else |
---|
1934 | printf ("854"); |
---|
1935 | if (IEEE) |
---|
1936 | printf (".\n"); |
---|
1937 | else { |
---|
1938 | printf (",\nexcept for possibly Double Rounding"); |
---|
1939 | printf (" during Gradual Underflow.\n"); |
---|
1940 | } |
---|
1941 | } |
---|
1942 | printf ("The arithmetic diagnosed appears to be Excellent!\n"); |
---|
1943 | } |
---|
1944 | } |
---|
1945 | |
---|
1946 | if (fpecount) |
---|
1947 | printf ("\nA total of %d floating point exceptions were registered.\n", |
---|
1948 | fpecount); |
---|
1949 | printf ("END OF TEST.\n"); |
---|
1950 | return 0; |
---|
1951 | } |
---|
1952 | |
---|
1953 | FLOAT |
---|
1954 | Sign (X) |
---|
1955 | FLOAT X; |
---|
1956 | { |
---|
1957 | return X >= 0. ? 1.0 : -1.0; |
---|
1958 | } |
---|
1959 | |
---|
1960 | void |
---|
1961 | Pause () |
---|
1962 | { |
---|
1963 | #ifndef BATCHMODE |
---|
1964 | char ch[8]; |
---|
1965 | |
---|
1966 | printf ("\nTo continue, press RETURN"); |
---|
1967 | fflush (stdout); |
---|
1968 | read (KEYBOARD, ch, 8); |
---|
1969 | #endif /* !BATCHMODE */ |
---|
1970 | #ifndef CYGNUS |
---|
1971 | printf ("\nDiagnosis resumes after milestone Number %d", Milestone); |
---|
1972 | printf (" Page: %d\n\n", PageNo); |
---|
1973 | ++Milestone; |
---|
1974 | ++PageNo; |
---|
1975 | #endif /* !CYGNUS */ |
---|
1976 | } |
---|
1977 | |
---|
1978 | void |
---|
1979 | TstCond (K, Valid, T) |
---|
1980 | int K, Valid; |
---|
1981 | char *T; |
---|
1982 | { |
---|
1983 | #ifdef CYGNUS |
---|
1984 | printf ("TEST: %s\n", T); |
---|
1985 | #endif /* CYGNUS */ |
---|
1986 | if (!Valid) { |
---|
1987 | BadCond (K, T); |
---|
1988 | printf (".\n"); |
---|
1989 | } |
---|
1990 | #ifdef CYGNUS |
---|
1991 | printf ("PASS: %s\n", T); |
---|
1992 | #endif /* CYGNUS */ |
---|
1993 | } |
---|
1994 | |
---|
1995 | void |
---|
1996 | BadCond (K, T) |
---|
1997 | int K; |
---|
1998 | char *T; |
---|
1999 | { |
---|
2000 | static char *msg[] = |
---|
2001 | {"FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW"}; |
---|
2002 | |
---|
2003 | ErrCnt[K] = ErrCnt[K] + 1; |
---|
2004 | #ifndef CYGNUS |
---|
2005 | printf ("%s: %s", msg[K], T); |
---|
2006 | #else |
---|
2007 | printf ("ERROR: Severity: %s: %s", msg[K], T); |
---|
2008 | #endif /* CYGNUS */ |
---|
2009 | } |
---|
2010 | |
---|
2011 | /* |
---|
2012 | * Random computes |
---|
2013 | * X = (Random1 + Random9)^5 |
---|
2014 | * Random1 = X - FLOOR(X) + 0.000005 * X; |
---|
2015 | * and returns the new value of Random1 |
---|
2016 | */ |
---|
2017 | FLOAT |
---|
2018 | Random () |
---|
2019 | { |
---|
2020 | FLOAT X, Y; |
---|
2021 | |
---|
2022 | X = Random1 + Random9; |
---|
2023 | Y = X * X; |
---|
2024 | Y = Y * Y; |
---|
2025 | X = X * Y; |
---|
2026 | Y = X - FLOOR (X); |
---|
2027 | Random1 = Y + X * 0.000005; |
---|
2028 | return (Random1); |
---|
2029 | } |
---|
2030 | |
---|
2031 | void |
---|
2032 | SqXMinX (ErrKind) |
---|
2033 | int ErrKind; |
---|
2034 | { |
---|
2035 | FLOAT XA, XB; |
---|
2036 | |
---|
2037 | XB = X * BInvrse; |
---|
2038 | XA = X - XB; |
---|
2039 | SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp; |
---|
2040 | if (SqEr != Zero) { |
---|
2041 | if (SqEr < MinSqEr) |
---|
2042 | MinSqEr = SqEr; |
---|
2043 | if (SqEr > MaxSqEr) |
---|
2044 | MaxSqEr = SqEr; |
---|
2045 | J = J + 1.0; |
---|
2046 | BadCond (ErrKind, "\n"); |
---|
2047 | printf ("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr); |
---|
2048 | printf ("\tinstead of correct value 0 .\n"); |
---|
2049 | } |
---|
2050 | } |
---|
2051 | |
---|
2052 | void |
---|
2053 | NewD () |
---|
2054 | { |
---|
2055 | X = Z1 * Q; |
---|
2056 | X = FLOOR (Half - X / Radix) * Radix + X; |
---|
2057 | Q = (Q - X * Z) / Radix + X * X * (D / Radix); |
---|
2058 | Z = Z - Two * X * D; |
---|
2059 | if (Z <= Zero) { |
---|
2060 | Z = -Z; |
---|
2061 | Z1 = -Z1; |
---|
2062 | } |
---|
2063 | D = Radix * D; |
---|
2064 | } |
---|
2065 | |
---|
2066 | void |
---|
2067 | SR3750 () |
---|
2068 | { |
---|
2069 | if (!((X - Radix < Z2 - Radix) || (X - Z2 > WVar - Z2))) { |
---|
2070 | I = I + 1; |
---|
2071 | X2 = SQRT (X * D); |
---|
2072 | Y2 = (X2 - Z2) - (Y - Z2); |
---|
2073 | X2 = X8 / (Y - Half); |
---|
2074 | X2 = X2 - Half * X2 * X2; |
---|
2075 | SqEr = (Y2 + Half) + (Half - X2); |
---|
2076 | if (SqEr < MinSqEr) |
---|
2077 | MinSqEr = SqEr; |
---|
2078 | SqEr = Y2 - X2; |
---|
2079 | if (SqEr > MaxSqEr) |
---|
2080 | MaxSqEr = SqEr; |
---|
2081 | } |
---|
2082 | } |
---|
2083 | |
---|
2084 | void |
---|
2085 | IsYeqX () |
---|
2086 | { |
---|
2087 | if (Y != X) { |
---|
2088 | if (N <= 0) { |
---|
2089 | if (Z == Zero && Q <= Zero) |
---|
2090 | printf ("WARNING: computing\n"); |
---|
2091 | else |
---|
2092 | BadCond (Defect, "computing\n"); |
---|
2093 | printf ("\t(%.17e) ^ (%.17e)\n", Z, Q); |
---|
2094 | printf ("\tyielded %.17e;\n", Y); |
---|
2095 | printf ("\twhich compared unequal to correct %.17e ;\n", |
---|
2096 | X); |
---|
2097 | printf ("\t\tthey differ by %.17e .\n", Y - X); |
---|
2098 | } |
---|
2099 | N = N + 1; /* ... count discrepancies. */ |
---|
2100 | } |
---|
2101 | } |
---|
2102 | |
---|
2103 | void |
---|
2104 | SR3980 () |
---|
2105 | { |
---|
2106 | do { |
---|
2107 | Q = (FLOAT) I; |
---|
2108 | Y = POW (Z, Q); |
---|
2109 | IsYeqX (); |
---|
2110 | if (++I > M) |
---|
2111 | break; |
---|
2112 | X = Z * X; |
---|
2113 | } |
---|
2114 | while (X < WVar); |
---|
2115 | } |
---|
2116 | |
---|
2117 | void |
---|
2118 | PrintIfNPositive () |
---|
2119 | { |
---|
2120 | if (N > 0) |
---|
2121 | printf ("Similar discrepancies have occurred %d times.\n", N); |
---|
2122 | } |
---|
2123 | |
---|
2124 | void |
---|
2125 | TstPtUf () |
---|
2126 | { |
---|
2127 | N = 0; |
---|
2128 | if (Z != Zero) { |
---|
2129 | printf ("Since comparison denies Z = 0, evaluating "); |
---|
2130 | printf ("(Z + Z) / Z should be safe.\n"); |
---|
2131 | sigsave = _sigfpe; |
---|
2132 | if (setjmp (ovfl_buf)) |
---|
2133 | goto very_serious; |
---|
2134 | Q9 = (Z + Z) / Z; |
---|
2135 | printf ("What the machine gets for (Z + Z) / Z is %.17e .\n", |
---|
2136 | Q9); |
---|
2137 | if (FABS (Q9 - Two) < Radix * U2) { |
---|
2138 | printf ("This is O.K., provided Over/Underflow"); |
---|
2139 | printf (" has NOT just been signaled.\n"); |
---|
2140 | } else { |
---|
2141 | if ((Q9 < One) || (Q9 > Two)) { |
---|
2142 | very_serious: |
---|
2143 | N = 1; |
---|
2144 | ErrCnt[Serious] = ErrCnt[Serious] + 1; |
---|
2145 | printf ("This is a VERY SERIOUS DEFECT!\n"); |
---|
2146 | } else { |
---|
2147 | N = 1; |
---|
2148 | ErrCnt[Defect] = ErrCnt[Defect] + 1; |
---|
2149 | printf ("This is a DEFECT!\n"); |
---|
2150 | } |
---|
2151 | } |
---|
2152 | sigsave = 0; |
---|
2153 | V9 = Z * One; |
---|
2154 | Random1 = V9; |
---|
2155 | V9 = One * Z; |
---|
2156 | Random2 = V9; |
---|
2157 | V9 = Z / One; |
---|
2158 | if ((Z == Random1) && (Z == Random2) && (Z == V9)) { |
---|
2159 | if (N > 0) |
---|
2160 | Pause (); |
---|
2161 | } else { |
---|
2162 | N = 1; |
---|
2163 | BadCond (Defect, "What prints as Z = "); |
---|
2164 | printf ("%.17e\n\tcompares different from ", Z); |
---|
2165 | if (Z != Random1) |
---|
2166 | printf ("Z * 1 = %.17e ", Random1); |
---|
2167 | if (!((Z == Random2) |
---|
2168 | || (Random2 == Random1))) |
---|
2169 | printf ("1 * Z == %g\n", Random2); |
---|
2170 | if (!(Z == V9)) |
---|
2171 | printf ("Z / 1 = %.17e\n", V9); |
---|
2172 | if (Random2 != Random1) { |
---|
2173 | ErrCnt[Defect] = ErrCnt[Defect] + 1; |
---|
2174 | BadCond (Defect, "Multiplication does not commute!\n"); |
---|
2175 | printf ("\tComparison alleges that 1 * Z = %.17e\n", |
---|
2176 | Random2); |
---|
2177 | printf ("\tdiffers from Z * 1 = %.17e\n", Random1); |
---|
2178 | } |
---|
2179 | Pause (); |
---|
2180 | } |
---|
2181 | } |
---|
2182 | } |
---|
2183 | |
---|
2184 | void |
---|
2185 | notify (s) |
---|
2186 | char *s; |
---|
2187 | { |
---|
2188 | printf ("%s test appears to be inconsistent...\n", s); |
---|
2189 | printf (" PLEASE NOTIFY KARPINKSI!\n"); |
---|
2190 | } |
---|
2191 | |
---|
2192 | void |
---|
2193 | msglist (s) |
---|
2194 | char **s; |
---|
2195 | { |
---|
2196 | while (*s) |
---|
2197 | printf ("%s\n", *s++); |
---|
2198 | } |
---|
2199 | |
---|
2200 | void |
---|
2201 | Instructions () |
---|
2202 | { |
---|
2203 | static char *instr[] = |
---|
2204 | { |
---|
2205 | "Lest this program stop prematurely, i.e. before displaying\n", |
---|
2206 | " `END OF TEST',\n", |
---|
2207 | "try to persuade the computer NOT to terminate execution when an", |
---|
2208 | "error like Over/Underflow or Division by Zero occurs, but rather", |
---|
2209 | "to persevere with a surrogate value after, perhaps, displaying some", |
---|
2210 | "warning. If persuasion avails naught, don't despair but run this", |
---|
2211 | "program anyway to see how many milestones it passes, and then", |
---|
2212 | "amend it to make further progress.\n", |
---|
2213 | "Answer questions with Y, y, N or n (unless otherwise indicated).\n", |
---|
2214 | 0}; |
---|
2215 | |
---|
2216 | msglist (instr); |
---|
2217 | } |
---|
2218 | |
---|
2219 | void |
---|
2220 | Heading () |
---|
2221 | { |
---|
2222 | static char *head[] = |
---|
2223 | { |
---|
2224 | "Users are invited to help debug and augment this program so it will", |
---|
2225 | "cope with unanticipated and newly uncovered arithmetic pathologies.\n", |
---|
2226 | "Please send suggestions and interesting results to", |
---|
2227 | "\tRichard Karpinski", |
---|
2228 | "\tComputer Center U-76", |
---|
2229 | "\tUniversity of California", |
---|
2230 | "\tSan Francisco, CA 94143-0704, USA\n", |
---|
2231 | "In doing so, please include the following information:", |
---|
2232 | #ifdef SINGLE_PRECISION |
---|
2233 | "\tPrecision:\tsingle;", |
---|
2234 | #else /* !SINGLE_PRECISION */ |
---|
2235 | "\tPrecision:\tdouble;", |
---|
2236 | #endif /* SINGLE_PRECISION */ |
---|
2237 | "\tVersion:\t10 February 1989;", |
---|
2238 | "\tComputer:\n", |
---|
2239 | "\tCompiler:\n", |
---|
2240 | "\tOptimization level:\n", |
---|
2241 | "\tOther relevant compiler options:", |
---|
2242 | 0}; |
---|
2243 | |
---|
2244 | msglist (head); |
---|
2245 | } |
---|
2246 | |
---|
2247 | void |
---|
2248 | Characteristics () |
---|
2249 | { |
---|
2250 | static char *chars[] = |
---|
2251 | { |
---|
2252 | "Running this program should reveal these characteristics:", |
---|
2253 | " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...", |
---|
2254 | " Precision = number of significant digits carried.", |
---|
2255 | " U2 = Radix/Radix^Precision = One Ulp", |
---|
2256 | "\t(OneUlpnit in the Last Place) of 1.000xxx .", |
---|
2257 | " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .", |
---|
2258 | " Adequacy of guard digits for Mult., Div. and Subt.", |
---|
2259 | " Whether arithmetic is chopped, correctly rounded, or something else", |
---|
2260 | "\tfor Mult., Div., Add/Subt. and Sqrt.", |
---|
2261 | " Whether a Sticky Bit used correctly for rounding.", |
---|
2262 | " UnderflowThreshold = an underflow threshold.", |
---|
2263 | " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.", |
---|
2264 | " V = an overflow threshold, roughly.", |
---|
2265 | " V0 tells, roughly, whether Infinity is represented.", |
---|
2266 | " Comparisions are checked for consistency with subtraction", |
---|
2267 | "\tand for contamination with pseudo-zeros.", |
---|
2268 | " Sqrt is tested. Y^X is not tested.", |
---|
2269 | " Extra-precise subexpressions are revealed but NOT YET tested.", |
---|
2270 | " Decimal-Binary conversion is NOT YET tested for accuracy.", |
---|
2271 | 0}; |
---|
2272 | |
---|
2273 | msglist (chars); |
---|
2274 | } |
---|
2275 | |
---|
2276 | void |
---|
2277 | History () |
---|
2278 | { /* History */ |
---|
2279 | /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner, |
---|
2280 | with further massaging by David M. Gay. */ |
---|
2281 | |
---|
2282 | static char *hist[] = |
---|
2283 | { |
---|
2284 | "The program attempts to discriminate among", |
---|
2285 | " FLAWs, like lack of a sticky bit,", |
---|
2286 | " Serious DEFECTs, like lack of a guard digit, and", |
---|
2287 | " FAILUREs, like 2+2 == 5 .", |
---|
2288 | "Failures may confound subsequent diagnoses.\n", |
---|
2289 | "The diagnostic capabilities of this program go beyond an earlier", |
---|
2290 | "program called `MACHAR', which can be found at the end of the", |
---|
2291 | "book `Software Manual for the Elementary Functions' (1980) by", |
---|
2292 | "W. J. Cody and W. Waite. Although both programs try to discover", |
---|
2293 | "the Radix, Precision and range (over/underflow thresholds)", |
---|
2294 | "of the arithmetic, this program tries to cope with a wider variety", |
---|
2295 | "of pathologies, and to say how well the arithmetic is implemented.", |
---|
2296 | "\nThe program is based upon a conventional radix representation for", |
---|
2297 | "floating-point numbers, but also allows logarithmic encoding", |
---|
2298 | "as used by certain early WANG machines.\n", |
---|
2299 | "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;", |
---|
2300 | "see source comments for more history.", |
---|
2301 | 0}; |
---|
2302 | |
---|
2303 | msglist (hist); |
---|
2304 | } |
---|