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