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