1 /* A C version of Kahan's Floating Point Test "Paranoia"
3 Thos Sumner, UCSF, Feb. 1985
4 David Gay, BTL, Jan. 1986
6 This is a rewrite from the Pascal version by
8 B. A. Wichmann, 18 Jan. 1985
10 (and does NOT exhibit good C programming style).
12 Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
13 compile with -DKR_headers or insert
15 at the beginning if you have an old-style C compiler.
17 (C) Apr 19 1983 in BASIC version by:
18 Professor W. M. Kahan,
20 Electrical Engineering & Computer Science Dept.
21 University of California
22 Berkeley, California 94720
25 converted to Pascal by:
27 National Physical Laboratory
34 David M. Gay and Thos Sumner
35 AT&T Bell Labs Computer Center, Rm. U-76
36 600 Mountain Avenue University of California
37 Murray Hill, NJ 07974 San Francisco, CA 94143
40 with simultaneous corrections to the Pascal source (reflected
41 in the Pascal source available over netlib).
42 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
44 Reports of results on various systems from all the versions
45 of Paranoia are being collected by Richard Karpinski at the
46 same address as Thos Sumner. This includes sample outputs,
47 bug reports, and criticisms.
49 You may copy this program freely if you acknowledge its source.
50 Comments on the Pascal version to NPL, please.
53 The C version catches signals from floating-point exceptions.
54 If signal(SIGFPE,...) is unavailable in your environment, you may
55 #define NOSIGNAL to comment out the invocations of signal.
57 This source file is too big for some C compilers, but may be split
58 into pieces. Comments containing "SPLIT" suggest convenient places
59 for this splitting. At the end of these comments is an "ed script"
60 (for the UNIX(tm) editor ed) that will do this splitting.
62 By #defining Single when you compile this source, you may obtain
63 a single-precision C version of Paranoia.
66 The following is from the introductory commentary from Wichmann's work:
68 The BASIC program of Kahan is written in Microsoft BASIC using many
69 facilities which have no exact analogy in Pascal. The Pascal
70 version below cannot therefore be exactly the same. Rather than be
71 a minimal transcription of the BASIC program, the Pascal coding
72 follows the conventional style of block-structured languages. Hence
73 the Pascal version could be useful in producing versions in other
76 Rather than use identifiers of minimal length (which therefore have
77 little mnemonic significance), the Pascal version uses meaningful
78 identifiers as follows [Note: A few changes have been made for C]:
81 BASIC C BASIC C BASIC C
84 A1 AInverse J0 NoErrors T
85 B Radix [Failure] T0 Underflow
86 B1 BInverse J1 NoErrors T2 ThirtyTwo
87 B2 RadixD2 [SeriousDefect] T5 OneAndHalf
88 B9 BMinusU2 J2 NoErrors T7 TwentySeven
89 C [Defect] T8 TwoForty
90 C1 CInverse J3 NoErrors U OneUlp
91 D [Flaw] U0 UnderflowThreshold
102 F1 MinusOne O5 Five X9 Random1
108 G2 GDiv Q9 Z0 PseudoZero
111 H1 HInverse R2 RDiv Z9
118 All the variables in BASIC are true variables and in consequence,
119 the program is more difficult to follow since the "constants" must
120 be determined (the glossary is very helpful). The Pascal version
121 uses Real constants, but checks are added to ensure that the values
122 are correctly converted by the compiler.
124 The major textual change to the Pascal version apart from the
125 identifiersis that named procedures are used, inserting parameters
126 wherehelpful. New procedures are also introduced. The
127 correspondence is as follows:
134 170- 250 Instructions
136 480- 670 Characteristics
140 4040-4080 DoesYequalX
141 4090-4110 PrintIfNPositive
142 4640-4850 TestPartialUnderflow
144 =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
146 Below is an "ed script" that splits para.c into 10 files
147 of the form part[1-8].c, subs.c, and msgs.c, plus a header
148 file, paranoia.h, that these files require.
207 /Computed constants/,$d
209 1,$s/^FLOAT/extern &/
212 /^Guard/,/^Round/s/^/extern /
213 /^jmp_buf/s/^/extern /
214 /^Sig_type/s/^/extern /
216 extern void sigfpe(INT);/
230 #define FABS(x) (float)fabs((double)(x))
231 #define FLOOR(x) (float)floor((double)(x))
232 #define LOG(x) (float)log((double)(x))
233 #define POW(x,y) (float)pow((double)(x),(double)(y))
234 #define SQRT(x) (float)sqrt((double)(x))
237 #define FABS(x) fabs(x)
238 #define FLOOR(x) floor(x)
239 #define LOG(x) log(x)
240 #define POW(x,y) pow(x,y)
241 #define SQRT(x) sqrt(x)
246 #define VOID /* void */
247 #define INT /* int */
248 #define FP /* FLOAT */
249 #define CHARP /* char * */
250 #define CHARPP /* char ** */
251 extern double fabs(), floor(), log(), pow(), sqrt();
253 typedef void (*Sig_type)();
254 FLOAT Sign(), Random();
255 extern void BadCond();
256 extern void SqXMinX();
257 extern void TstCond();
258 extern void notify();
265 #define CHARPP char **
273 extern double fabs(double), floor(double), log(double);
274 extern double pow(double,double), sqrt(double);
275 extern void exit(INT);
280 typedef void (*Sig_type)(int);
281 FLOAT Sign(FLOAT), Random(void);
282 extern void BadCond(int, char*);
283 extern void SqXMinX(int);
284 extern void TstCond(int, int, char*);
285 extern void notify(char*);
286 extern int read(int, char*, int);
289 extern void Characteristics(VOID);
290 extern void Heading(VOID);
291 extern void History(VOID);
292 extern void Instructions(VOID);
293 extern void IsYeqX(VOID);
294 extern void NewD(VOID);
295 extern void Pause(VOID);
296 extern void PrintIfNPositive(VOID);
297 extern void SR3750(VOID);
298 extern void SR3980(VOID);
299 extern void TstPtUf(VOID);
305 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
307 /*Small floating point constants.*/
317 FLOAT TwentySeven = 27.0;
318 FLOAT ThirtyTwo = 32.0;
319 FLOAT TwoForty = 240.0;
320 FLOAT MinusOne = -1.0;
321 FLOAT OneAndHalf = 1.5;
322 /*Integer constants*/
323 int NoTrials = 20; /*Number of tests for commutativity. */
327 /* Definitions for declared types
329 Rounding == (Chopped, Rounded, Other);
330 Message == packed array [1..40] of char;
331 Class == (Flaw, Defect, Serious, Failure);
342 typedef int Guard, Rounding, Class;
343 typedef char Message;
345 /* Declarations of Variables */
351 FLOAT E0, E1, Exp2, E3, MinSqEr;
352 FLOAT SqEr, MaxSqEr, E9;
362 FLOAT T, Underflow, S;
363 FLOAT OneUlp, UfThold, U1, U2;
366 FLOAT X, X1, X2, X8, Random1;
367 FLOAT Y, Y1, Y2, Random2;
368 FLOAT Z, PseudoZero, Z1, Z2, Z9;
374 Guard GMult, GDiv, GAddSub;
375 Rounding RMult, RDiv, RAddSub, RSqrt;
376 int Break, Done, NotMonot, Monot, Anomaly, IEEE,
378 /* Computed constants. */
379 /*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
380 /*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
382 /* floating point exception receiver */
387 printf("\n* * * FLOATING-POINT ERROR %d * * *\n", x);
391 signal(SIGFPE, sigsave);
394 longjmp(ovfl_buf, 1);
401 /* First two assignments use integer right-hand sides. */
409 Nine = Three * Three;
410 TwentySeven = Nine * Three;
411 ThirtyTwo = Four * Eight;
412 TwoForty = Four * Five * Three * Four;
415 OneAndHalf = One + Half;
421 /*=============================================*/
423 /*=============================================*/
425 signal(SIGFPE, sigfpe);
435 /*=============================================*/
437 /*=============================================*/
438 printf("Program is now RUNNING tests on small integers:\n");
440 TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
441 && (One > Zero) && (One + One == Two),
442 "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
445 ErrCnt[Failure] = ErrCnt[Failure] + 1;
446 printf("Comparison alleges that -0.0 is Non-zero!\n");
451 TstCond (Failure, (Three == Two + One) && (Four == Three + One)
452 && (Four + Two * (- Two) == Zero)
453 && (Four - Three - One == Zero),
454 "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
455 TstCond (Failure, (MinusOne == (0 - One))
456 && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
457 && (MinusOne + FABS(One) == Zero)
458 && (MinusOne + MinusOne * MinusOne == Zero),
459 "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
460 TstCond (Failure, Half + MinusOne + Half == Zero,
461 "1/2 + (-1) + 1/2 != 0");
462 /*=============================================*/
465 extern void part2(VOID), part3(VOID), part4(VOID),
466 part5(VOID), part6(VOID), part7(VOID);
478 #include "paranoia.h"
482 /*=============================================*/
483 TstCond (Failure, (Nine == Three * Three)
484 && (TwentySeven == Nine * Three) && (Eight == Four + Four)
485 && (ThirtyTwo == Eight * Four)
486 && (ThirtyTwo - TwentySeven - Four - One == Zero),
487 "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
488 TstCond (Failure, (Five == Four + One) &&
489 (TwoForty == Four * Five * Three * Four)
490 && (TwoForty / Three - Four * Four * Five == Zero)
491 && ( TwoForty / Four - Five * Three * Four == Zero)
492 && ( TwoForty / Five - Four * Three * Four == Zero),
493 "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
494 if (ErrCnt[Failure] == 0) {
495 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
498 printf("Searching for Radix and Precision.\n");
505 } while (MinusOne + FABS(Y) < Zero);
506 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
513 } while ( Radix == Zero);
514 if (Radix < Two) Radix = One;
515 printf("Radix = %f .\n", Radix);
519 Precision = Precision + One;
522 } while ((Y - W) == One);
524 /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
528 printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
529 printf("Recalculating radix and precision\n ");
544 /*... now X = (unknown no.) ulps of 1+...*/
547 Y = Half * U2 + ThirtyTwo * U2 * U2;
550 } while ( ! ((U2 <= X) || (X <= Zero)));
552 /*... now U2 == 1 ulp of 1 + ... */
560 /*... now X == (unknown no.) ulps of 1 -... */
563 Y = Half * U1 + ThirtyTwo * U1 * U1;
568 } while ( ! ((U1 <= X) || (X <= Zero)));
569 /*... now U1 == 1 ulp of 1 - ... */
570 if (U1 == E1) printf("confirms closest relative separation U1 .\n");
571 else printf("gets better closest relative separation U1 = %.7e .\n", U1);
573 F9 = (Half - U1) + Half;
574 Radix = FLOOR(0.01 + U2 / U1);
575 if (Radix == E0) printf("Radix confirmed.\n");
576 else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
577 TstCond (Defect, Radix <= Eight + Eight,
578 "Radix is too big: roundoff problems");
579 TstCond (Flaw, (Radix == Two) || (Radix == 10)
580 || (Radix == One), "Radix is not as good as 2 or 10");
581 /*=============================================*/
583 /*=============================================*/
584 TstCond (Failure, F9 - Half < Half,
585 "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
590 TstCond (Failure, (X != One)
591 || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
594 /*=============================================*/
596 /*=============================================*/
597 /*... BMinusU2 = nextafter(Radix, 0) */
598 BMinusU2 = Radix - One;
599 BMinusU2 = (BMinusU2 - U2) + One;
600 /* Purify Integers */
602 X = - TwoForty * LOG(U1) / LOG(Radix);
604 if (FABS(X - Y) * Four < One) X = Y;
605 Precision = X / TwoForty;
606 Y = FLOOR(Half + Precision);
607 if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
609 if ((Precision != FLOOR(Precision)) || (Radix == One)) {
610 printf("Precision cannot be characterized by an Integer number\n");
611 printf("of significant digits but, by itself, this is a minor flaw.\n");
614 printf("logarithmic encoding has precision characterized solely by U1.\n");
615 else printf("The number of significant digits of the Radix is %f .\n",
617 TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
618 "Precision worse than 5 decimal figures ");
619 /*=============================================*/
621 /*=============================================*/
622 /* Test for extra-precise subepressions */
623 X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
626 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
627 } while ( ! ((Z2 <= X) || (X <= Zero)));
628 X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
631 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
632 + One / Two)) + One / Two;
633 } while ( ! ((Z1 <= Z) || (Z <= Zero)));
637 Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
639 } while ( ! ((Y1 <= Y) || (Y <= Zero)));
641 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
642 } while ( ! ((X1 <= X) || (X <= Zero)));
643 if ((X1 != Y1) || (X1 != Z1)) {
644 BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
645 printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
646 printf("are symptoms of inconsistencies introduced\n");
647 printf("by extra-precise evaluation of arithmetic subexpressions.\n");
648 notify("Possibly some part of this");
649 if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
650 "That feature is not tested further by this program.\n") ;
653 if ((Z1 != U1) || (Z2 != U2)) {
654 if ((Z1 >= U1) || (Z2 >= U2)) {
655 BadCond(Failure, "");
657 printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
658 printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
661 if ((Z1 <= Zero) || (Z2 <= Zero)) {
662 printf("Because of unusual Radix = %f", Radix);
663 printf(", or exact rational arithmetic a result\n");
664 printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
665 notify("of an\nextra-precision");
667 if (Z1 != Z2 || Z1 > Zero) {
672 printf("Some subexpressions appear to be calculated extra\n");
673 printf("precisely with about %g extra B-digits, i.e.\n",
675 printf("roughly %g extra significant decimals.\n",
678 printf("That feature is not tested further by this program.\n");
683 /*=============================================*/
686 #include "paranoia.h"
690 /*=============================================*/
692 X = W / (Radix * Radix);
697 TstCond (Failure, X == U2,
698 "Subtraction is not normalized X=Y,X+Z != Y+Z!");
700 "Subtraction appears to be normalized, as it should be.");
702 printf("\nChecking for guard digit in *, /, and -.\n");
715 X = X * (Radix - One);
716 T = T * (Radix - One);
717 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
720 TstCond (Serious, False,
721 "* lacks a Guard Digit, so 1*X != X");
725 Y = FABS((X + Z) - X * X) - U2;
727 Z = FABS((X - U2) - X * X) - U1;
728 TstCond (Failure, (Y <= Zero)
729 && (Z <= Zero), "* gets too many final digits wrong.\n");
737 T = Nine / TwentySeven;
739 TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
740 "Division lacks a Guard Digit, so error can exceed 1 ulp\n\
741 or 1/3 and 3/9 and 9/27 may disagree");
748 if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
751 TstCond (Serious, False,
752 "Division lacks a Guard Digit, so X/1 != X");
754 X = One / (One + U2);
756 TstCond (Serious, Y < Zero,
757 "Computed value of 1/1.000..1 >= 1");
759 Y = One + Radix * U2;
763 StickyBit = T / Radix;
766 TstCond (Failure, X == Zero && Y == Zero,
767 "* and/or / gets too many last digits wrong");
772 Z = Radix - BMinusU2;
774 if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
777 TstCond (Serious, False,
778 "- lacks Guard Digit, so cancellation is obscured");
780 if (F9 != One && F9 - One >= Zero) {
781 BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
782 printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
783 printf(" such precautions against division by zero as\n");
784 printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
786 if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
787 " *, /, and - appear to have guard digits, as they should.\n");
788 /*=============================================*/
790 /*=============================================*/
792 printf("Checking rounding on multiply, divide and add/subtract.\n");
796 RadixD2 = Radix / Two;
803 AInvrse = AInvrse / A1;
804 } while ( ! (FLOOR(AInvrse) != AInvrse));
805 Done = (X == One) || (A1 > Three);
806 if (! Done) A1 = Nine + One;
808 if (X == One) A1 = Radix;
815 TstCond (Failure, Z == Half,
816 "X * (1/X) differs from 1");
833 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
834 X = (OneAndHalf + U2) * Y2;
835 Y = OneAndHalf - U2 - U2;
836 Z = OneAndHalf + U2 + U2;
837 T = (OneAndHalf - U2) * Y1;
842 Y = (U2 - Y) + StickyBit;
843 Z = S - (Z + U2 + U2);
844 StickyBit = (Y2 + U2) * Y1;
846 StickyBit = StickyBit - Y2;
848 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
849 && ( StickyBit == Zero) && (Y1 == Half)) {
851 printf("Multiplication appears to round correctly.\n");
853 else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
854 && (T < Zero) && (StickyBit + U2 == Zero)
857 printf("Multiplication appears to chop.\n");
859 else printf("* is neither chopped nor correctly rounded.\n");
860 if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
862 else printf("* is neither chopped nor correctly rounded.\n");
863 /*=============================================*/
865 /*=============================================*/
868 Z = OneAndHalf + U2 + U2;
870 T = OneAndHalf - U2 - U2;
876 Z = Z - (OneAndHalf + U2);
877 T = (U2 - OneAndHalf) + T;
878 if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
892 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
893 && (Y2 == Zero) && (Y2 == Zero)
894 && (Y1 - Half == F9 - Half )) {
896 printf("Division appears to round correctly.\n");
897 if (GDiv == No) notify("Division");
899 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
900 && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
902 printf("Division appears to chop.\n");
905 if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
906 BInvrse = One / Radix;
907 TstCond (Failure, (BInvrse * Radix - Half == Half),
908 "Radix * ( 1 / Radix ) differs from 1");
909 /*=============================================*/
912 #include "paranoia.h"
916 /*=============================================*/
917 TstCond (Failure, ((F9 + U1) - Half == Half)
918 && ((BMinusU2 + U2 ) - One == Radix - One),
919 "Incomplete carry-propagation in Addition");
921 Y = One + U2 * (One - U2);
925 if ((X == Zero) && (Y == Zero)) {
927 printf("Add/Subtract appears to be chopped.\n");
929 if (GAddSub == Yes) {
930 X = (Half + U2) * U2;
931 Y = (Half - U2) * U2;
936 if ((X == Zero) && (Y == Zero)) {
937 X = (Half + U2) * U1;
938 Y = (Half - U2) * U1;
943 if ((X == Zero) && (Y == Zero)) {
945 printf("Addition/Subtraction appears to round correctly.\n");
946 if (GAddSub == No) notify("Add/Subtract");
948 else printf("Addition/Subtraction neither rounds nor chops.\n");
950 else printf("Addition/Subtraction neither rounds nor chops.\n");
952 else printf("Addition/Subtraction neither rounds nor chops.\n");
954 X = One + Half * (One + Half);
955 Y = (One + U2) * Half;
959 if (StickyBit != Zero) {
961 BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
964 if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
965 && (RMult == Rounded) && (RDiv == Rounded)
966 && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
967 printf("Checking for sticky bit.\n");
968 X = (Half + U1) * U2;
972 if ((Z - One <= Zero) && (T - One >= U2)) {
975 if ((Z - T >= U2) && (Y - T == Zero)) {
976 X = (Half + U1) * U1;
980 if ((Z - One == Zero) && (T - F9 == Zero)) {
981 Z = (Half - U1) * U1;
984 if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
985 Z = (One + U2) * OneAndHalf;
986 T = (OneAndHalf + U2) - Z + U2;
987 X = One + Half / Radix;
988 Y = One + Radix * U2;
990 if (T == Zero && X + Radix * U2 - Z == Zero) {
994 if ((Y - One == Zero)) StickyBit = S;
1003 if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
1004 else printf("Sticky bit used incorrectly or not at all.\n");
1005 TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
1006 RMult == Other || RDiv == Other || RAddSub == Other),
1007 "lack(s) of guard digits or failure(s) to correctly round or chop\n\
1008 (noted above) count as one flaw in the final tally below");
1009 /*=============================================*/
1011 /*=============================================*/
1013 printf("Does Multiplication commute? ");
1014 printf("Testing on %d random pairs.\n", NoTrials);
1015 Random9 = SQRT(3.0);
1025 } while ( ! ((I > NoTrials) || (Z9 != Zero)));
1026 if (I == NoTrials) {
1027 Random1 = One + Half / Three;
1028 Random2 = (U2 + U1) + One;
1029 Z = Random1 * Random2;
1030 Y = Random2 * Random1;
1031 Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
1032 Three) * ((U2 + U1) + One);
1034 if (! ((I == NoTrials) || (Z9 == Zero)))
1035 BadCond(Defect, "X * Y == Y * X trial fails.\n");
1036 else printf(" No failures found in %d integer pairs.\n", NoTrials);
1037 /*=============================================*/
1039 /*=============================================*/
1040 printf("\nRunning test of square root(x).\n");
1041 TstCond (Failure, (Zero == SQRT(Zero))
1042 && (- Zero == SQRT(- Zero))
1043 && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
1051 OneUlp = BInvrse * U1;
1056 if (J != Zero) Pause();
1057 printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
1061 if ((Radix != One)) do {
1064 } while ( ! ((Y - X >= NoTrials)));
1067 while (I <= NoTrials) {
1070 if (J > Zero) break;
1073 printf("Test for sqrt monotonicity.\n");
1077 Z = Radix + Radix * U2;
1080 while ( ! (NotMonot || Monot)) {
1085 if ((X > Q) || (Q > Z)) NotMonot = True;
1087 Q = FLOOR(Q + Half);
1088 if (!(I > 0 || Radix == Q * Q)) Monot = True;
1090 if (I > 1) Monot = True;
1104 if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
1106 BadCond(Defect, "");
1107 printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
1109 /*=============================================*/
1112 #include "paranoia.h"
1116 /*=============================================*/
1117 MinSqEr = MinSqEr + Half;
1118 MaxSqEr = MaxSqEr - Half;
1119 Y = (SQRT(One + U2) - One) / U2;
1120 SqEr = (Y - One) + U2 / Eight;
1121 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1122 SqEr = Y + U2 / Eight;
1123 if (SqEr < MinSqEr) MinSqEr = SqEr;
1124 Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
1125 SqEr = Y + U1 / Eight;
1126 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1127 SqEr = (Y + One) + U1 / Eight;
1128 if (SqEr < MinSqEr) MinSqEr = SqEr;
1131 for( Indx = 1; Indx <= 3; ++Indx) {
1132 Y = SQRT((X + U1 + X) + F9);
1133 Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
1134 Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
1135 SqEr = (Y + Half) + Z;
1136 if (SqEr < MinSqEr) MinSqEr = SqEr;
1137 SqEr = (Y - Half) + Z;
1138 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1139 if (((Indx == 1) || (Indx == 3)))
1140 X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
1146 /*=============================================*/
1148 /*=============================================*/
1151 RSqrt = Other; /* ~dgh */
1153 printf("Testing whether sqrt is rounded or chopped.\n");
1154 D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
1155 /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
1158 if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
1173 X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
1176 } while ( ! (X1 <= Zero));
1185 if (Y2 >= FourD) Y2 = Y2 - FourD;
1186 } while ( ! (Y >= D));
1188 Q = (X8 + Z * Z) / FourD;
1190 if (Q != FLOOR(Q)) Anomaly = True;
1195 X = X - FLOOR(X / Radix) * Radix;
1200 } while ( ! (Break || (Z1 <= Zero)));
1201 if ((Z1 <= Zero) && (! Break)) Anomaly = True;
1203 if (Z1 > RadixD2) Z1 = Z1 - Radix;
1206 } while ( ! (U2 * D >= F9));
1207 if (D * Radix - D != W - D) Anomaly = True;
1211 Y = D + (One + Z) * Half;
1214 Y = D + (One - Z) * Half + D;
1219 if (D - Z2 != W - Z2) Anomaly = True;
1221 Y = (D - Z2) + (Z2 + (One - Z) * Half);
1222 X = (D - Z2) + (Z2 - Z + Q);
1224 Y = (One + Z) * Half;
1227 if (I == 0) Anomaly = True;
1233 if ((I == 0) || Anomaly) {
1234 BadCond(Failure, "Anomalous arithmetic with Integer < ");
1235 printf("Radix^Precision = %.7e\n", W);
1236 printf(" fails test whether sqrt rounds or chops.\n");
1241 if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
1243 printf("Square root appears to be correctly rounded.\n");
1246 if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
1247 || (MinSqEr + Radix < Half)) SqRWrng = True;
1250 printf("Square root appears to be chopped.\n");
1255 printf("Square root is neither chopped nor correctly rounded.\n");
1256 printf("Observed errors run from %.7e ", MinSqEr - Half);
1257 printf("to %.7e ulps.\n", Half + MaxSqEr);
1258 TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
1259 "sqrt gets too many last digits wrong");
1261 /*=============================================*/
1263 /*=============================================*/
1265 printf("Testing powers Z^i for small Integers Z and i.\n");
1267 /* ... test powers of zero. */
1279 if (Z == MinusOne) Break = True;
1282 /* .. if(-1)^N is invalid, replace MinusOne by One. */
1290 M = (int)FLOOR(Two * LOG(W) / LOG(A1));
1296 if (Z == AInvrse) Break = True;
1298 } while ( ! (Break));
1299 /*=============================================*/
1301 /*=============================================*/
1302 /* Powers of Radix have been tested, */
1303 /* next try a few primes */
1312 } while ( Three * FLOOR(Z / Three) == Z );
1313 } while ( Z < Eight * Three );
1315 printf("Errors like this may invalidate financial calculations\n");
1316 printf("\tinvolving interest rates.\n");
1320 if (N == 0) printf("... no discrepancies found.\n");
1323 /*=============================================*/
1326 #include "paranoia.h"
1330 /*=============================================*/
1331 printf("Seeking Underflow thresholds UfThold and E0.\n");
1333 if (Precision != FLOOR(Precision)) {
1339 } while ( X > Zero);
1343 /* ... D is power of 1/Radix < 1. */
1348 } while ((Y > Z) && (Z + Z > Z));
1355 } while ((Y > Z) && (Z + Z > Z));
1356 if (Radix < Two) HInvrse = Two;
1357 else HInvrse = Radix;
1359 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1363 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1368 } while ((E0 > Z) && (Z + Z > Z));
1380 BadCond(Failure, "multiplication gets too many last digits wrong.\n");
1389 PseudoZero = Underflow * H;
1393 Underflow = PseudoZero;
1394 if (E1 + E1 <= E1) {
1395 Y2 = Underflow * HInvrse;
1398 if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
1400 PseudoZero = PseudoZero * H;
1401 } while ((Underflow > PseudoZero)
1402 && (PseudoZero + PseudoZero > PseudoZero));
1404 /* Comment line 4530 .. 4560 */
1405 if (PseudoZero != Zero) {
1408 /* ... Test PseudoZero for "phoney- zero" violates */
1409 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1411 if (PseudoZero <= Zero) {
1412 BadCond(Failure, "Positive expressions can underflow to an\n");
1413 printf("allegedly negative value\n");
1414 printf("PseudoZero that prints out as: %g .\n", PseudoZero);
1417 printf("But -PseudoZero, which should be\n");
1418 printf("positive, isn't; it prints out as %g .\n", X);
1422 BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
1423 printf("value PseudoZero that prints out as %g .\n", PseudoZero);
1427 /*=============================================*/
1429 /*=============================================*/
1430 if (CInvrse * Y > CInvrse * Y1) {
1434 if (! ((E1 == Zero) || (E1 == E0))) {
1435 BadCond(Defect, "");
1437 printf("Products underflow at a higher");
1438 printf(" threshold than differences.\n");
1439 if (PseudoZero == Zero)
1443 printf("Difference underflows at a higher");
1444 printf(" threshold than products.\n");
1447 printf("Smallest strictly positive number found is E0 = %g .\n", E0);
1451 if (N == 1) Underflow = Y;
1453 if (E1 == Zero) I = 3;
1454 if (UfThold == Zero) I = I - 2;
1458 UfThold = Underflow;
1459 if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
1461 BadCond(Failure, "Either accuracy deteriorates as numbers\n");
1462 printf("approach a threshold = %.17e\n", UfThold);;
1463 printf(" coming down from %.17e\n", C);
1464 printf(" or else multiplication gets too many last digits wrong.\n");
1470 BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
1471 printf("Q == Y while denying that |Q - Y| == 0; these values\n");
1472 printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
1473 printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
1482 if ((Q == UfThold) && (E1 == E0)
1483 && (FABS( UfThold - E1 / E9) <= E1)) {
1485 printf("Underflow is gradual; it incurs Absolute Error =\n");
1486 printf("(roundoff in UfThold) < E0.\n");
1488 Y = Y * (OneAndHalf + U2);
1489 X = CInvrse * (One + U2);
1497 if (setjmp(ovfl_buf)) {
1498 printf("Underflow / UfThold failed!\n");
1501 else R = SQRT(Underflow / UfThold);
1505 X = Z * (One + R * H * (One + H));
1509 X = Z * (One + H * H * (One + H));
1511 if (! ((X == Z) || (X - Z != Zero))) {
1513 printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
1515 printf("yet X - Z yields %.17e .\n", Z9);
1516 printf(" Should this NOT signal Underflow, ");
1517 printf("this is a SERIOUS DEFECT\nthat causes ");
1518 printf("confusion when innocent statements like\n");;
1519 printf(" if (X == Z) ... else");
1520 printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
1521 printf("encounter Division by Zero although actually\n");
1523 if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
1524 else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
1528 printf("The Underflow threshold is %.17e, %s\n", UfThold,
1530 printf("calculation may suffer larger Relative error than ");
1531 printf("merely roundoff.\n");
1535 if (Y2 <= UfThold) {
1537 BadCond(Defect, "");
1541 BadCond(Serious, "");
1544 printf("Range is too narrow; U1^%d Underflows.\n", I);
1546 /*=============================================*/
1549 #include "paranoia.h"
1553 /*=============================================*/
1554 Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
1556 printf("Since underflow occurs below the threshold\n");
1557 printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
1558 printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n",
1560 printf("actually calculating yields:");
1561 if (setjmp(ovfl_buf)) {
1563 BadCond(Serious, "trap on underflow.\n");
1567 V9 = POW(HInvrse, Y2);
1569 printf(" %.17e .\n", V9);
1570 if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
1571 BadCond(Serious, "this is not between 0 and underflow\n");
1572 printf(" threshold = %.17e .\n", UfThold);
1574 else if (! (V9 > UfThold * (One + E9)))
1575 printf("This computed value is O.K.\n");
1577 BadCond(Defect, "this is not between 0 and underflow\n");
1578 printf(" threshold = %.17e .\n", UfThold);
1581 /*=============================================*/
1583 /*=============================================*/
1585 /* ...calculate Exp2 == exp(2) == 7.389056099... */
1599 Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
1604 printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
1608 Z = (X + One) / (Z - (One - BInvrse));
1609 Q = POW(X, Z) - Exp2;
1610 if (FABS(Q) > TwoForty * U2) {
1612 V9 = (X - BInvrse) - (One - BInvrse);
1613 BadCond(Defect, "Calculated");
1614 printf(" %.17e for\n", POW(X,Z));
1615 printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
1616 printf("\tdiffers from correct value by %.17e .\n", Q);
1617 printf("\tThis much error may spoil financial\n");
1618 printf("\tcalculations involving tiny interest rates.\n");
1622 Z = (Y - X) * Two + Y;
1625 Z = One + (X - F9)*(X - F9);
1626 if (Z > One && I < NoTrials) I++;
1630 printf("Accuracy seems adequate.\n");
1642 /*=============================================*/
1644 /*=============================================*/
1645 printf("Testing powers Z^Q at four nearly extreme values.\n");
1648 Q = FLOOR(Half - LOG(C) / LOG(A1));
1658 if (Z < One) Break = True;
1660 } while ( ! (Break));
1662 if (N == 0) printf(" ... no discrepancies found.\n");
1665 /*=============================================*/
1667 /*=============================================*/
1669 printf("Searching for Overflow threshold:\n");
1670 printf("This may generate an error.\n");
1674 if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
1684 printf("Can `Z = -Y' overflow?\n");
1685 printf("Trying it on Y = %.17e .\n", Y);
1688 if (V - Y == V + V0) printf("Seems O.K.\n");
1691 BadCond(Flaw, "-(-Y) differs from Y.\n");
1694 BadCond(Serious, "");
1695 printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
1698 Y = V * (HInvrse * U2 - HInvrse);
1699 Z = Y + ((One - HInvrse) * U2) * V;
1702 if (V0 - V < V0) V = V0;
1705 V = Y * (HInvrse * U2 - HInvrse);
1706 V = V + ((One - HInvrse) * U2) * Y;
1708 printf("Overflow threshold is V = %.17e .\n", V);
1709 if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
1710 else printf("There is no saturation value because \
1711 the system traps on overflow.\n");
1713 printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
1715 printf(" nor for V / 1 = %.17e .\n", V9);
1716 printf("Any overflow signal separating this * from the one\n");
1717 printf("above is a DEFECT.\n");
1718 /*=============================================*/
1720 /*=============================================*/
1721 if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
1722 BadCond(Failure, "Comparisons involving ");
1723 printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
1726 /*=============================================*/
1728 /*=============================================*/
1730 for(Indx = 1; Indx <= 3; ++Indx) {
1732 case 1: Z = UfThold; break;
1733 case 2: Z = E0; break;
1734 case 3: Z = PseudoZero; break;
1739 if (Y / (One - Radix * E9) < Z
1740 || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
1741 if (V9 > U1) BadCond(Serious, "");
1742 else BadCond(Defect, "");
1743 printf("Comparison alleges that what prints as Z = %.17e\n", Z);
1744 printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
1748 /*=============================================*/
1750 /*=============================================*/
1751 for(Indx = 1; Indx <= 2; ++Indx) {
1752 if (Indx == 1) Z = V;
1755 X = (One - Radix * E9) * V9;
1757 if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
1759 if (X < W) BadCond(Serious, "");
1760 else BadCond(Defect, "");
1761 printf("Comparison alleges that Z = %17e\n", Z);
1762 printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
1765 /*=============================================*/
1768 #include "paranoia.h"
1772 /*=============================================*/
1776 if (X*Y < One || X > Y) {
1777 if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
1778 else BadCond(Flaw, "");
1780 printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
1781 X, "is too far from 1.\n");
1783 /*=============================================*/
1785 /*=============================================*/
1786 for (Indx = 1; Indx <= 5; ++Indx) {
1789 case 2: X = One + U2; break;
1790 case 3: X = V; break;
1791 case 4: X = UfThold; break;
1796 if (setjmp(ovfl_buf))
1797 printf(" X / X traps when X = %g\n", X);
1799 V9 = (Y / X - Half) - Half;
1800 if (V9 == Zero) continue;
1801 if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
1802 else BadCond(Serious, "");
1803 printf(" X / X differs from 1 when X = %.17e\n", X);
1804 printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
1808 /*=============================================*/
1810 /*=============================================*/
1813 printf("What message and/or values does Division by Zero produce?\n") ;
1815 printf("This can interupt your program. You can ");
1816 printf("skip this part if you wish.\n");
1817 printf("Do you wish to compute 1 / 0? ");
1819 read (KEYBOARD, ch, 8);
1820 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1823 printf(" Trying to compute 1 / 0 produces ...");
1824 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
1828 else printf("O.K.\n");
1829 printf("\nDo you wish to compute 0 / 0? ");
1831 read (KEYBOARD, ch, 80);
1832 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1835 printf("\n Trying to compute 0 / 0 produces ...");
1836 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
1840 else printf("O.K.\n");
1842 /*=============================================*/
1844 /*=============================================*/
1848 static char *msg[] = {
1849 "FAILUREs encountered =",
1850 "SERIOUS DEFECTs discovered =",
1851 "DEFECTs discovered =",
1852 "FLAWs discovered =" };
1854 for(i = 0; i < 4; i++) if (ErrCnt[i])
1855 printf("The number of %-29s %d.\n",
1859 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
1860 + ErrCnt[Flaw]) > 0) {
1861 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
1862 Defect] == 0) && (ErrCnt[Flaw] > 0)) {
1863 printf("The arithmetic diagnosed seems ");
1864 printf("Satisfactory though flawed.\n");
1866 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
1867 && ( ErrCnt[Defect] > 0)) {
1868 printf("The arithmetic diagnosed may be Acceptable\n");
1869 printf("despite inconvenient Defects.\n");
1871 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
1872 printf("The arithmetic diagnosed has ");
1873 printf("unacceptable Serious Defects.\n");
1875 if (ErrCnt[Failure] > 0) {
1876 printf("Potentially fatal FAILURE may have spoiled this");
1877 printf(" program's subsequent diagnoses.\n");
1881 printf("No failures, defects nor flaws have been discovered.\n");
1882 if (! ((RMult == Rounded) && (RDiv == Rounded)
1883 && (RAddSub == Rounded) && (RSqrt == Rounded)))
1884 printf("The arithmetic diagnosed seems Satisfactory.\n");
1886 if (StickyBit >= One &&
1887 (Radix - Two) * (Radix - Nine - One) == Zero) {
1888 printf("Rounding appears to conform to ");
1889 printf("the proposed IEEE standard P");
1890 if ((Radix == Two) &&
1891 ((Precision - Four * Three * Two) *
1892 ( Precision - TwentySeven -
1893 TwentySeven + One) == Zero))
1896 if (IEEE) printf(".\n");
1898 printf(",\nexcept for possibly Double Rounding");
1899 printf(" during Gradual Underflow.\n");
1902 printf("The arithmetic diagnosed appears to be Excellent!\n");
1906 printf("\nA total of %d floating point exceptions were registered.\n",
1908 printf("END OF TEST.\n");
1913 #include "paranoia.h"
1921 { return X >= 0. ? 1.0 : -1.0; }
1929 printf("\nTo continue, press RETURN");
1931 read(KEYBOARD, ch, 8);
1933 printf("\nDiagnosis resumes after milestone Number %d", Milestone);
1934 printf(" Page: %d\n\n", PageNo);
1940 TstCond (INT K, INT Valid, CHARP T)
1945 { if (! Valid) { BadCond(K,T); printf(".\n"); } }
1948 BadCond(INT K, CHARP T)
1954 static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
1956 ErrCnt [K] = ErrCnt [K] + 1;
1957 printf("%s: %s", msg[K], T);
1964 X = (Random1 + Random9)^5
1965 Random1 = X - FLOOR(X) + 0.000005 * X;
1966 and returns the new value of Random1
1971 X = Random1 + Random9;
1976 Random1 = Y + X * 0.000005;
1981 SqXMinX (INT ErrKind)
1990 SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
1992 if (SqEr < MinSqEr) MinSqEr = SqEr;
1993 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1995 BadCond(ErrKind, "\n");
1996 printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
1997 printf("\tinstead of correct value 0 .\n");
2005 X = FLOOR(Half - X / Radix) * Radix + X;
2006 Q = (Q - X * Z) / Radix + X * X * (D / Radix);
2007 Z = Z - Two * X * D;
2018 if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
2021 Y2 = (X2 - Z2) - (Y - Z2);
2022 X2 = X8 / (Y - Half);
2023 X2 = X2 - Half * X2 * X2;
2024 SqEr = (Y2 + Half) + (Half - X2);
2025 if (SqEr < MinSqEr) MinSqEr = SqEr;
2027 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
2036 if (Z == Zero && Q <= Zero)
2037 printf("WARNING: computing\n");
2038 else BadCond(Defect, "computing\n");
2039 printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
2040 printf("\tyielded %.17e;\n", Y);
2041 printf("\twhich compared unequal to correct %.17e ;\n",
2043 printf("\t\tthey differ by %.17e .\n", Y - X);
2045 N = N + 1; /* ... count discrepancies. */
2062 PrintIfNPositive(VOID)
2064 if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
2072 printf("Since comparison denies Z = 0, evaluating ");
2073 printf("(Z + Z) / Z should be safe.\n");
2075 if (setjmp(ovfl_buf)) goto very_serious;
2077 printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
2079 if (FABS(Q9 - Two) < Radix * U2) {
2080 printf("This is O.K., provided Over/Underflow");
2081 printf(" has NOT just been signaled.\n");
2084 if ((Q9 < One) || (Q9 > Two)) {
2087 ErrCnt [Serious] = ErrCnt [Serious] + 1;
2088 printf("This is a VERY SERIOUS DEFECT!\n");
2092 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2093 printf("This is a DEFECT!\n");
2102 if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
2107 BadCond(Defect, "What prints as Z = ");
2108 printf("%.17e\n\tcompares different from ", Z);
2109 if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
2110 if (! ((Z == Random2)
2111 || (Random2 == Random1)))
2112 printf("1 * Z == %g\n", Random2);
2113 if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
2114 if (Random2 != Random1) {
2115 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2116 BadCond(Defect, "Multiplication does not commute!\n");
2117 printf("\tComparison alleges that 1 * Z = %.17e\n",
2119 printf("\tdiffers from Z * 1 = %.17e\n", Random1);
2132 printf("%s test appears to be inconsistent...\n", s);
2133 printf(" PLEASE NOTIFY KARPINKSI!\n");
2137 #include "paranoia.h"
2145 { while(*s) printf("%s\n", *s++); }
2150 static char *instr[] = {
2151 "Lest this program stop prematurely, i.e. before displaying\n",
2152 " `END OF TEST',\n",
2153 "try to persuade the computer NOT to terminate execution when an",
2154 "error like Over/Underflow or Division by Zero occurs, but rather",
2155 "to persevere with a surrogate value after, perhaps, displaying some",
2156 "warning. If persuasion avails naught, don't despair but run this",
2157 "program anyway to see how many milestones it passes, and then",
2158 "amend it to make further progress.\n",
2159 "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
2168 static char *head[] = {
2169 "Users are invited to help debug and augment this program so it will",
2170 "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
2171 "Please send suggestions and interesting results to",
2172 "\tRichard Karpinski",
2173 "\tComputer Center U-76",
2174 "\tUniversity of California",
2175 "\tSan Francisco, CA 94143-0704, USA\n",
2176 "In doing so, please include the following information:",
2178 "\tPrecision:\tsingle;",
2180 "\tPrecision:\tdouble;",
2182 "\tVersion:\t10 February 1989;",
2185 "\tOptimization level:\n",
2186 "\tOther relevant compiler options:",
2193 Characteristics(VOID)
2195 static char *chars[] = {
2196 "Running this program should reveal these characteristics:",
2197 " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
2198 " Precision = number of significant digits carried.",
2199 " U2 = Radix/Radix^Precision = One Ulp",
2200 "\t(OneUlpnit in the Last Place) of 1.000xxx .",
2201 " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
2202 " Adequacy of guard digits for Mult., Div. and Subt.",
2203 " Whether arithmetic is chopped, correctly rounded, or something else",
2204 "\tfor Mult., Div., Add/Subt. and Sqrt.",
2205 " Whether a Sticky Bit used correctly for rounding.",
2206 " UnderflowThreshold = an underflow threshold.",
2207 " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
2208 " V = an overflow threshold, roughly.",
2209 " V0 tells, roughly, whether Infinity is represented.",
2210 " Comparisions are checked for consistency with subtraction",
2211 "\tand for contamination with pseudo-zeros.",
2212 " Sqrt is tested. Y^X is not tested.",
2213 " Extra-precise subexpressions are revealed but NOT YET tested.",
2214 " Decimal-Binary conversion is NOT YET tested for accuracy.",
2223 /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
2224 with further massaging by David M. Gay. */
2226 static char *hist[] = {
2227 "The program attempts to discriminate among",
2228 " FLAWs, like lack of a sticky bit,",
2229 " Serious DEFECTs, like lack of a guard digit, and",
2230 " FAILUREs, like 2+2 == 5 .",
2231 "Failures may confound subsequent diagnoses.\n",
2232 "The diagnostic capabilities of this program go beyond an earlier",
2233 "program called `MACHAR', which can be found at the end of the",
2234 "book `Software Manual for the Elementary Functions' (1980) by",
2235 "W. J. Cody and W. Waite. Although both programs try to discover",
2236 "the Radix, Precision and range (over/underflow thresholds)",
2237 "of the arithmetic, this program tries to cope with a wider variety",
2238 "of pathologies, and to say how well the arithmetic is implemented.",
2239 "\nThe program is based upon a conventional radix representation for",
2240 "floating-point numbers, but also allows logarithmic encoding",
2241 "as used by certain early WANG machines.\n",
2242 "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
2243 "see source comments for more history.",