2 /* A C version of Kahan's Floating Point Test "Paranoia"
4 Thos Sumner, UCSF, Feb. 1985
5 David Gay, BTL, Jan. 1986
7 This is a rewrite from the Pascal version by
9 B. A. Wichmann, 18 Jan. 1985
11 (and does NOT exhibit good C programming style).
13 Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
14 compile with -DKR_headers or insert
16 at the beginning if you have an old-style C compiler.
18 (C) Apr 19 1983 in BASIC version by:
19 Professor W. M. Kahan,
21 Electrical Engineering & Computer Science Dept.
22 University of California
23 Berkeley, California 94720
26 converted to Pascal by:
28 National Physical Laboratory
35 David M. Gay and Thos Sumner
36 AT&T Bell Labs Computer Center, Rm. U-76
37 600 Mountain Avenue University of California
38 Murray Hill, NJ 07974 San Francisco, CA 94143
41 with simultaneous corrections to the Pascal source (reflected
42 in the Pascal source available over netlib).
43 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
45 Reports of results on various systems from all the versions
46 of Paranoia are being collected by Richard Karpinski at the
47 same address as Thos Sumner. This includes sample outputs,
48 bug reports, and criticisms.
50 You may copy this program freely if you acknowledge its source.
51 Comments on the Pascal version to NPL, please.
54 The C version catches signals from floating-point exceptions.
55 If signal(SIGFPE,...) is unavailable in your environment, you may
56 #define NOSIGNAL to comment out the invocations of signal.
58 This source file is too big for some C compilers, but may be split
59 into pieces. Comments containing "SPLIT" suggest convenient places
60 for this splitting. At the end of these comments is an "ed script"
61 (for the UNIX(tm) editor ed) that will do this splitting.
63 By #defining Single when you compile this source, you may obtain
64 a single-precision C version of Paranoia.
67 The following is from the introductory commentary from Wichmann's work:
69 The BASIC program of Kahan is written in Microsoft BASIC using many
70 facilities which have no exact analogy in Pascal. The Pascal
71 version below cannot therefore be exactly the same. Rather than be
72 a minimal transcription of the BASIC program, the Pascal coding
73 follows the conventional style of block-structured languages. Hence
74 the Pascal version could be useful in producing versions in other
77 Rather than use identifiers of minimal length (which therefore have
78 little mnemonic significance), the Pascal version uses meaningful
79 identifiers as follows [Note: A few changes have been made for C]:
82 BASIC C BASIC C BASIC C
85 A1 AInverse J0 NoErrors T
86 B Radix [Failure] T0 Underflow
87 B1 BInverse J1 NoErrors T2 ThirtyTwo
88 B2 RadixD2 [SeriousDefect] T5 OneAndHalf
89 B9 BMinusU2 J2 NoErrors T7 TwentySeven
90 C [Defect] T8 TwoForty
91 C1 CInverse J3 NoErrors U OneUlp
92 D [Flaw] U0 UnderflowThreshold
103 F1 MinusOne O5 Five X9 Random1
109 G2 GDiv Q9 Z0 PseudoZero
112 H1 HInverse R2 RDiv Z9
119 All the variables in BASIC are true variables and in consequence,
120 the program is more difficult to follow since the "constants" must
121 be determined (the glossary is very helpful). The Pascal version
122 uses Real constants, but checks are added to ensure that the values
123 are correctly converted by the compiler.
125 The major textual change to the Pascal version apart from the
126 identifiersis that named procedures are used, inserting parameters
127 wherehelpful. New procedures are also introduced. The
128 correspondence is as follows:
135 170- 250 Instructions
137 480- 670 Characteristics
141 4040-4080 DoesYequalX
142 4090-4110 PrintIfNPositive
143 4640-4850 TestPartialUnderflow
145 =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
147 Below is an "ed script" that splits para.c into 10 files
148 of the form part[1-8].c, subs.c, and msgs.c, plus a header
149 file, paranoia.h, that these files require.
208 /Computed constants/,$d
210 1,$s/^FLOAT/extern &/
213 /^Guard/,/^Round/s/^/extern /
214 /^jmp_buf/s/^/extern /
215 /^Sig_type/s/^/extern /
217 extern void sigfpe(INT);/
229 #if #cpu(r4640) || #cpu(r4650)
235 #define FABS(x) fabsf((x))
236 #define FLOOR(x) (float)floor((double)(x))
237 #define LOG(x) (float)log((double)(x))
238 #define POW(x,y) (float)pow((double)(x),(double)(y))
239 #define SQRT(x) sqrtf((x))
242 #define FABS(x) fabs(x)
243 #define FLOOR(x) floor(x)
244 #define LOG(x) log(x)
245 #define POW(x,y) pow(x,y)
246 #define SQRT(x) sqrt(x)
251 #define VOID /* void */
252 #define INT /* int */
253 #define FP /* FLOAT */
254 #define CHARP /* char * */
255 #define CHARPP /* char ** */
256 extern double fabs(), floor(), log(), pow(), sqrt();
258 typedef void (*Sig_type)();
259 FLOAT Sign(), Random();
260 extern void BadCond();
261 extern void SqXMinX();
262 extern void TstCond();
263 extern void notify();
270 #define CHARPP char **
278 extern double fabs(double), floor(double), log(double);
279 extern double pow(double,double), sqrt(double);
280 extern void exit(INT);
285 typedef void (*Sig_type)(int);
286 FLOAT Sign(FLOAT), Random(void);
287 extern void BadCond(int, char*);
288 extern void SqXMinX(int);
289 extern void TstCond(int, int, char*);
290 extern void notify(char*);
291 extern int read(int, char*, int);
294 extern void Characteristics(VOID);
295 extern void Heading(VOID);
296 extern void History(VOID);
297 extern void Instructions(VOID);
298 extern void IsYeqX(VOID);
299 extern void NewD(VOID);
300 extern void Pause(VOID);
301 extern void PrintIfNPositive(VOID);
302 extern void SR3750(VOID);
303 extern void SR3980(VOID);
304 extern void TstPtUf(VOID);
310 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
312 /*Small floating point constants.*/
322 FLOAT TwentySeven = 27.0;
323 FLOAT ThirtyTwo = 32.0;
324 FLOAT TwoForty = 240.0;
325 FLOAT MinusOne = -1.0;
326 FLOAT OneAndHalf = 1.5;
327 /*Integer constants*/
328 int NoTrials = 20; /*Number of tests for commutativity. */
332 /* Definitions for declared types
334 Rounding == (Chopped, Rounded, Other);
335 Message == packed array [1..40] of char;
336 Class == (Flaw, Defect, Serious, Failure);
347 typedef int Guard, Rounding, Class;
348 typedef char Message;
350 /* Declarations of Variables */
356 FLOAT E0, E1, Exp2, E3, MinSqEr;
357 FLOAT SqEr, MaxSqEr, E9;
367 FLOAT T, Underflow, S;
368 FLOAT OneUlp, UfThold, U1, U2;
371 FLOAT X, X1, X2, X8, Random1;
372 FLOAT Y, Y1, Y2, Random2;
373 FLOAT Z, PseudoZero, Z1, Z2, Z9;
379 Guard GMult, GDiv, GAddSub;
380 Rounding RMult, RDiv, RAddSub, RSqrt;
381 int Break, Done, NotMonot, Monot, Anomaly, IEEE,
383 /* Computed constants. */
384 /*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
385 /*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
387 /* floating point exception receiver */
392 printf("\n* * * FLOATING-POINT ERROR %d * * *\n", x);
396 signal(SIGFPE, sigsave);
399 longjmp(ovfl_buf, 1);
407 /* First two assignments use integer right-hand sides. */
415 Nine = Three * Three;
416 TwentySeven = Nine * Three;
417 ThirtyTwo = Four * Eight;
418 TwoForty = Four * Five * Three * Four;
421 OneAndHalf = One + Half;
427 /*=============================================*/
429 /*=============================================*/
431 signal(SIGFPE, sigfpe);
441 /*=============================================*/
443 /*=============================================*/
444 printf("Program is now RUNNING tests on small integers:\n");
446 TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
447 && (One > Zero) && (One + One == Two),
448 "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
450 if (Z != (FLOAT)0.0) {
451 ErrCnt[Failure] = ErrCnt[Failure] + 1;
452 printf("Comparison alleges that -0.0 is Non-zero!\n");
457 TstCond (Failure, (Three == Two + One) && (Four == Three + One)
458 && (Four + Two * (- Two) == Zero)
459 && (Four - Three - One == Zero),
460 "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
461 TstCond (Failure, (MinusOne == (0 - One))
462 && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
463 && (MinusOne + FABS(One) == Zero)
464 && (MinusOne + MinusOne * MinusOne == Zero),
465 "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
466 TstCond (Failure, Half + MinusOne + Half == Zero,
467 "1/2 + (-1) + 1/2 != 0");
468 /*=============================================*/
471 extern void part2(VOID), part3(VOID), part4(VOID),
472 part5(VOID), part6(VOID), part7(VOID);
484 #include "paranoia.h"
488 /*=============================================*/
489 TstCond (Failure, (Nine == Three * Three)
490 && (TwentySeven == Nine * Three) && (Eight == Four + Four)
491 && (ThirtyTwo == Eight * Four)
492 && (ThirtyTwo - TwentySeven - Four - One == Zero),
493 "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
494 TstCond (Failure, (Five == Four + One) &&
495 (TwoForty == Four * Five * Three * Four)
496 && (TwoForty / Three - Four * Four * Five == Zero)
497 && ( TwoForty / Four - Five * Three * Four == Zero)
498 && ( TwoForty / Five - Four * Three * Four == Zero),
499 "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
500 if (ErrCnt[Failure] == 0) {
501 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
504 printf("Searching for Radix and Precision.\n");
511 } while (MinusOne + FABS(Y) < Zero);
512 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
519 } while ( Radix == Zero);
520 if (Radix < Two) Radix = One;
521 printf("Radix = %f .\n", Radix);
525 Precision = Precision + One;
528 } while ((Y - W) == One);
530 /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
534 printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
535 printf("Recalculating radix and precision\n ");
550 /*... now X = (unknown no.) ulps of 1+...*/
553 Y = Half * U2 + ThirtyTwo * U2 * U2;
556 } while ( ! ((U2 <= X) || (X <= Zero)));
558 /*... now U2 == 1 ulp of 1 + ... */
566 /*... now X == (unknown no.) ulps of 1 -... */
569 Y = Half * U1 + ThirtyTwo * U1 * U1;
574 } while ( ! ((U1 <= X) || (X <= Zero)));
575 /*... now U1 == 1 ulp of 1 - ... */
576 if (U1 == E1) printf("confirms closest relative separation U1 .\n");
577 else printf("gets better closest relative separation U1 = %.7e .\n", U1);
579 F9 = (Half - U1) + Half;
580 Radix = FLOOR((FLOAT)0.01 + U2 / U1);
581 if (Radix == E0) printf("Radix confirmed.\n");
582 else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
583 TstCond (Defect, Radix <= Eight + Eight,
584 "Radix is too big: roundoff problems");
585 TstCond (Flaw, (Radix == Two) || (Radix == 10)
586 || (Radix == One), "Radix is not as good as 2 or 10");
587 /*=============================================*/
589 /*=============================================*/
590 TstCond (Failure, F9 - Half < Half,
591 "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
596 TstCond (Failure, (X != One)
597 || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
600 /*=============================================*/
602 /*=============================================*/
603 /*... BMinusU2 = nextafter(Radix, 0) */
604 BMinusU2 = Radix - One;
605 BMinusU2 = (BMinusU2 - U2) + One;
606 /* Purify Integers */
608 X = - TwoForty * LOG(U1) / LOG(Radix);
610 if (FABS(X - Y) * Four < One) X = Y;
611 Precision = X / TwoForty;
612 Y = FLOOR(Half + Precision);
613 if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
615 if ((Precision != FLOOR(Precision)) || (Radix == One)) {
616 printf("Precision cannot be characterized by an Integer number\n");
617 printf("of significant digits but, by itself, this is a minor flaw.\n");
620 printf("logarithmic encoding has precision characterized solely by U1.\n");
621 else printf("The number of significant digits of the Radix is %f .\n",
623 TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
624 "Precision worse than 5 decimal figures ");
625 /*=============================================*/
627 /*=============================================*/
628 /* Test for extra-precise subepressions */
629 X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
632 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
633 } while ( ! ((Z2 <= X) || (X <= Zero)));
634 X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
637 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
638 + One / Two)) + One / Two;
639 } while ( ! ((Z1 <= Z) || (Z <= Zero)));
643 Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
645 } while ( ! ((Y1 <= Y) || (Y <= Zero)));
647 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
648 } while ( ! ((X1 <= X) || (X <= Zero)));
649 if ((X1 != Y1) || (X1 != Z1)) {
650 BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
651 printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
652 printf("are symptoms of inconsistencies introduced\n");
653 printf("by extra-precise evaluation of arithmetic subexpressions.\n");
654 notify("Possibly some part of this");
655 if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
656 "That feature is not tested further by this program.\n") ;
659 if ((Z1 != U1) || (Z2 != U2)) {
660 if ((Z1 >= U1) || (Z2 >= U2)) {
661 BadCond(Failure, "");
663 printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
664 printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
667 if ((Z1 <= Zero) || (Z2 <= Zero)) {
668 printf("Because of unusual Radix = %f", Radix);
669 printf(", or exact rational arithmetic a result\n");
670 printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
671 notify("of an\nextra-precision");
673 if (Z1 != Z2 || Z1 > Zero) {
678 printf("Some subexpressions appear to be calculated extra\n");
679 printf("precisely with about %g extra B-digits, i.e.\n",
681 printf("roughly %g extra significant decimals.\n",
684 printf("That feature is not tested further by this program.\n");
689 /*=============================================*/
692 #include "paranoia.h"
696 /*=============================================*/
698 X = W / (Radix * Radix);
703 TstCond (Failure, X == U2,
704 "Subtraction is not normalized X=Y,X+Z != Y+Z!");
706 "Subtraction appears to be normalized, as it should be.");
708 printf("\nChecking for guard digit in *, /, and -.\n");
721 X = X * (Radix - One);
722 T = T * (Radix - One);
723 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
726 TstCond (Serious, False,
727 "* lacks a Guard Digit, so 1*X != X");
731 Y = FABS((X + Z) - X * X) - U2;
733 Z = FABS((X - U2) - X * X) - U1;
734 TstCond (Failure, (Y <= Zero)
735 && (Z <= Zero), "* gets too many final digits wrong.\n");
743 T = Nine / TwentySeven;
745 TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
746 "Division lacks a Guard Digit, so error can exceed 1 ulp\n\
747 or 1/3 and 3/9 and 9/27 may disagree");
754 if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
757 TstCond (Serious, False,
758 "Division lacks a Guard Digit, so X/1 != X");
760 X = One / (One + U2);
762 TstCond (Serious, Y < Zero,
763 "Computed value of 1/1.000..1 >= 1");
765 Y = One + Radix * U2;
769 StickyBit = T / Radix;
772 TstCond (Failure, X == Zero && Y == Zero,
773 "* and/or / gets too many last digits wrong");
778 Z = Radix - BMinusU2;
780 if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
783 TstCond (Serious, False,
784 "- lacks Guard Digit, so cancellation is obscured");
786 if (F9 != One && F9 - One >= Zero) {
787 BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
788 printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
789 printf(" such precautions against division by zero as\n");
790 printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
792 if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
793 " *, /, and - appear to have guard digits, as they should.\n");
794 /*=============================================*/
796 /*=============================================*/
798 printf("Checking rounding on multiply, divide and add/subtract.\n");
802 RadixD2 = Radix / Two;
809 AInvrse = AInvrse / A1;
810 } while ( ! (FLOOR(AInvrse) != AInvrse));
811 Done = (X == One) || (A1 > Three);
812 if (! Done) A1 = Nine + One;
814 if (X == One) A1 = Radix;
821 TstCond (Failure, Z == Half,
822 "X * (1/X) differs from 1");
839 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
840 X = (OneAndHalf + U2) * Y2;
841 Y = OneAndHalf - U2 - U2;
842 Z = OneAndHalf + U2 + U2;
843 T = (OneAndHalf - U2) * Y1;
848 Y = (U2 - Y) + StickyBit;
849 Z = S - (Z + U2 + U2);
850 StickyBit = (Y2 + U2) * Y1;
852 StickyBit = StickyBit - Y2;
854 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
855 && ( StickyBit == Zero) && (Y1 == Half)) {
857 printf("Multiplication appears to round correctly.\n");
859 else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
860 && (T < Zero) && (StickyBit + U2 == Zero)
863 printf("Multiplication appears to chop.\n");
865 else printf("* is neither chopped nor correctly rounded.\n");
866 if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
868 else printf("* is neither chopped nor correctly rounded.\n");
869 /*=============================================*/
871 /*=============================================*/
874 Z = OneAndHalf + U2 + U2;
876 T = OneAndHalf - U2 - U2;
882 Z = Z - (OneAndHalf + U2);
883 T = (U2 - OneAndHalf) + T;
884 if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
898 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
899 && (Y2 == Zero) && (Y2 == Zero)
900 && (Y1 - Half == F9 - Half )) {
902 printf("Division appears to round correctly.\n");
903 if (GDiv == No) notify("Division");
905 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
906 && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
908 printf("Division appears to chop.\n");
911 if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
912 BInvrse = One / Radix;
913 TstCond (Failure, (BInvrse * Radix - Half == Half),
914 "Radix * ( 1 / Radix ) differs from 1");
915 /*=============================================*/
918 #include "paranoia.h"
922 /*=============================================*/
923 TstCond (Failure, ((F9 + U1) - Half == Half)
924 && ((BMinusU2 + U2 ) - One == Radix - One),
925 "Incomplete carry-propagation in Addition");
927 Y = One + U2 * (One - U2);
931 if ((X == Zero) && (Y == Zero)) {
933 printf("Add/Subtract appears to be chopped.\n");
935 if (GAddSub == Yes) {
936 X = (Half + U2) * U2;
937 Y = (Half - U2) * U2;
942 if ((X == Zero) && (Y == Zero)) {
943 X = (Half + U2) * U1;
944 Y = (Half - U2) * U1;
949 if ((X == Zero) && (Y == Zero)) {
951 printf("Addition/Subtraction appears to round correctly.\n");
952 if (GAddSub == No) notify("Add/Subtract");
954 else printf("Addition/Subtraction neither rounds nor chops.\n");
956 else printf("Addition/Subtraction neither rounds nor chops.\n");
958 else printf("Addition/Subtraction neither rounds nor chops.\n");
960 X = One + Half * (One + Half);
961 Y = (One + U2) * Half;
965 if (StickyBit != Zero) {
967 BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
970 if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
971 && (RMult == Rounded) && (RDiv == Rounded)
972 && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
973 printf("Checking for sticky bit.\n");
974 X = (Half + U1) * U2;
978 if ((Z - One <= Zero) && (T - One >= U2)) {
981 if ((Z - T >= U2) && (Y - T == Zero)) {
982 X = (Half + U1) * U1;
986 if ((Z - One == Zero) && (T - F9 == Zero)) {
987 Z = (Half - U1) * U1;
990 if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
991 Z = (One + U2) * OneAndHalf;
992 T = (OneAndHalf + U2) - Z + U2;
993 X = One + Half / Radix;
994 Y = One + Radix * U2;
996 if (T == Zero && X + Radix * U2 - Z == Zero) {
1000 if ((Y - One == Zero)) StickyBit = S;
1009 if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
1010 else printf("Sticky bit used incorrectly or not at all.\n");
1011 TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
1012 RMult == Other || RDiv == Other || RAddSub == Other),
1013 "lack(s) of guard digits or failure(s) to correctly round or chop\n\
1014 (noted above) count as one flaw in the final tally below");
1015 /*=============================================*/
1017 /*=============================================*/
1019 printf("Does Multiplication commute? ");
1020 printf("Testing on %d random pairs.\n", NoTrials);
1021 Random9 = SQRT(3.0);
1031 } while ( ! ((I > NoTrials) || (Z9 != Zero)));
1032 if (I == NoTrials) {
1033 Random1 = One + Half / Three;
1034 Random2 = (U2 + U1) + One;
1035 Z = Random1 * Random2;
1036 Y = Random2 * Random1;
1037 Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
1038 Three) * ((U2 + U1) + One);
1040 if (! ((I == NoTrials) || (Z9 == Zero)))
1041 BadCond(Defect, "X * Y == Y * X trial fails.\n");
1042 else printf(" No failures found in %d integer pairs.\n", NoTrials);
1043 /*=============================================*/
1045 /*=============================================*/
1046 printf("\nRunning test of square root(x).\n");
1047 TstCond (Failure, (Zero == SQRT(Zero))
1048 && (- Zero == SQRT(- Zero))
1049 && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
1057 OneUlp = BInvrse * U1;
1062 if (J != Zero) Pause();
1063 printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
1067 if ((Radix != One)) do {
1070 } while ( ! ((Y - X >= NoTrials)));
1073 while (I <= NoTrials) {
1076 if (J > Zero) break;
1079 printf("Test for sqrt monotonicity.\n");
1083 Z = Radix + Radix * U2;
1086 while ( ! (NotMonot || Monot)) {
1091 if ((X > Q) || (Q > Z)) NotMonot = True;
1093 Q = FLOOR(Q + Half);
1094 if ((I > 0) || (Radix == Q * Q)) Monot = True;
1096 if (I > 1) Monot = True;
1110 if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
1112 BadCond(Defect, "");
1113 printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
1115 /*=============================================*/
1118 #include "paranoia.h"
1122 /*=============================================*/
1123 MinSqEr = MinSqEr + Half;
1124 MaxSqEr = MaxSqEr - Half;
1125 Y = (SQRT(One + U2) - One) / U2;
1126 SqEr = (Y - One) + U2 / Eight;
1127 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1128 SqEr = Y + U2 / Eight;
1129 if (SqEr < MinSqEr) MinSqEr = SqEr;
1130 Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
1131 SqEr = Y + U1 / Eight;
1132 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1133 SqEr = (Y + One) + U1 / Eight;
1134 if (SqEr < MinSqEr) MinSqEr = SqEr;
1137 for( Indx = 1; Indx <= 3; ++Indx) {
1138 Y = SQRT((X + U1 + X) + F9);
1139 Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
1140 Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
1141 SqEr = (Y + Half) + Z;
1142 if (SqEr < MinSqEr) MinSqEr = SqEr;
1143 SqEr = (Y - Half) + Z;
1144 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1145 if (((Indx == 1) || (Indx == 3)))
1146 X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
1152 /*=============================================*/
1154 /*=============================================*/
1157 RSqrt = Other; /* ~dgh */
1159 printf("Testing whether sqrt is rounded or chopped.\n");
1160 D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
1161 /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
1164 if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
1179 X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
1182 } while ( ! (X1 <= Zero));
1191 if (Y2 >= FourD) Y2 = Y2 - FourD;
1192 } while ( ! (Y >= D));
1194 Q = (X8 + Z * Z) / FourD;
1196 if (Q != FLOOR(Q)) Anomaly = True;
1201 X = X - FLOOR(X / Radix) * Radix;
1206 } while ( ! (Break || (Z1 <= Zero)));
1207 if ((Z1 <= Zero) && (! Break)) Anomaly = True;
1209 if (Z1 > RadixD2) Z1 = Z1 - Radix;
1212 } while ( ! (U2 * D >= F9));
1213 if (D * Radix - D != W - D) Anomaly = True;
1217 Y = D + (One + Z) * Half;
1220 Y = D + (One - Z) * Half + D;
1225 if (D - Z2 != W - Z2) Anomaly = True;
1227 Y = (D - Z2) + (Z2 + (One - Z) * Half);
1228 X = (D - Z2) + (Z2 - Z + Q);
1230 Y = (One + Z) * Half;
1233 if (I == 0) Anomaly = True;
1239 if ((I == 0) || Anomaly) {
1240 BadCond(Failure, "Anomalous arithmetic with Integer < ");
1241 printf("Radix^Precision = %.7e\n", W);
1242 printf(" fails test whether sqrt rounds or chops.\n");
1247 if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
1249 printf("Square root appears to be correctly rounded.\n");
1252 if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
1253 || (MinSqEr + Radix < Half)) SqRWrng = True;
1256 printf("Square root appears to be chopped.\n");
1261 printf("Square root is neither chopped nor correctly rounded.\n");
1262 printf("Observed errors run from %.7e ", MinSqEr - Half);
1263 printf("to %.7e ulps.\n", Half + MaxSqEr);
1264 TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
1265 "sqrt gets too many last digits wrong");
1267 /*=============================================*/
1269 /*=============================================*/
1271 printf("Testing powers Z^i for small Integers Z and i.\n");
1273 /* ... test powers of zero. */
1285 if (Z == MinusOne) Break = True;
1288 /* .. if(-1)^N is invalid, replace MinusOne by One. */
1296 M = (int)FLOOR(Two * LOG(W) / LOG(A1));
1302 if (Z == AInvrse) Break = True;
1304 } while ( ! (Break));
1305 /*=============================================*/
1307 /*=============================================*/
1308 /* Powers of Radix have been tested, */
1309 /* next try a few primes */
1318 } while ( Three * FLOOR(Z / Three) == Z );
1319 } while ( Z < Eight * Three );
1321 printf("Errors like this may invalidate financial calculations\n");
1322 printf("\tinvolving interest rates.\n");
1326 if (N == 0) printf("... no discrepancies found.\n");
1329 /*=============================================*/
1332 #include "paranoia.h"
1336 /*=============================================*/
1337 printf("Seeking Underflow thresholds UfThold and E0.\n");
1339 if (Precision != FLOOR(Precision)) {
1345 } while ( X > Zero);
1349 printf("Y: %20.20f Z: %20.20f\n", Y, Z);
1350 /* ... D is power of 1/Radix < 1. */
1355 } while ((Y > Z) && (Z + Z > Z));
1362 } while ((Y > Z) && (Z + Z > Z));
1363 if (Radix < Two) HInvrse = Two;
1364 else HInvrse = Radix;
1366 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1370 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1375 } while ((E0 > Z) && (Z + Z > Z));
1387 BadCond(Failure, "multiplication gets too many last digits wrong.\n");
1396 PseudoZero = Underflow * H;
1400 Underflow = PseudoZero;
1401 if (E1 + E1 <= E1) {
1402 Y2 = Underflow * HInvrse;
1405 if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
1407 PseudoZero = PseudoZero * H;
1408 } while ((Underflow > PseudoZero)
1409 && (PseudoZero + PseudoZero > PseudoZero));
1411 /* Comment line 4530 .. 4560 */
1412 if (PseudoZero != Zero) {
1415 /* ... Test PseudoZero for "phoney- zero" violates */
1416 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1418 if (PseudoZero <= Zero) {
1419 BadCond(Failure, "Positive expressions can underflow to an\n");
1420 printf("allegedly negative value\n");
1421 printf("PseudoZero that prints out as: %g .\n", PseudoZero);
1424 printf("But -PseudoZero, which should be\n");
1425 printf("positive, isn't; it prints out as %g .\n", X);
1429 BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
1430 printf("value PseudoZero that prints out as %g .\n", PseudoZero);
1434 /*=============================================*/
1436 /*=============================================*/
1437 if (CInvrse * Y > CInvrse * Y1) {
1441 if (! ((E1 == Zero) || (E1 == E0))) {
1442 BadCond(Defect, "");
1444 printf("Products underflow at a higher");
1445 printf(" threshold than differences.\n");
1446 if (PseudoZero == Zero)
1450 printf("Difference underflows at a higher");
1451 printf(" threshold than products.\n");
1454 printf("Smallest strictly positive number found is E0 = %g .\n", E0);
1458 if (N == 1) Underflow = Y;
1460 if (E1 == Zero) I = 3;
1461 if (UfThold == Zero) I = I - 2;
1465 UfThold = Underflow;
1466 if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
1468 BadCond(Failure, "Either accuracy deteriorates as numbers\n");
1469 printf("approach a threshold = %.17e\n", UfThold);;
1470 printf(" coming down from %.17e\n", C);
1471 printf(" or else multiplication gets too many last digits wrong.\n");
1477 BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
1478 printf("Q == Y while denying that |Q - Y| == 0; these values\n");
1479 printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
1480 printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
1489 if ((Q == UfThold) && (E1 == E0)
1490 && (FABS( UfThold - E1 / E9) <= E1)) {
1492 printf("Underflow is gradual; it incurs Absolute Error =\n");
1493 printf("(roundoff in UfThold) < E0.\n");
1495 Y = Y * (OneAndHalf + U2);
1496 X = CInvrse * (One + U2);
1504 if (setjmp(ovfl_buf)) {
1505 printf("Underflow / UfThold failed!\n");
1508 else R = SQRT(Underflow / UfThold);
1512 X = Z * (One + R * H * (One + H));
1516 X = Z * (One + H * H * (One + H));
1518 if (! ((X == Z) || (X - Z != Zero))) {
1520 printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
1522 printf("yet X - Z yields %.17e .\n", Z9);
1523 printf(" Should this NOT signal Underflow, ");
1524 printf("this is a SERIOUS DEFECT\nthat causes ");
1525 printf("confusion when innocent statements like\n");;
1526 printf(" if (X == Z) ... else");
1527 printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
1528 printf("encounter Division by Zero although actually\n");
1530 if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
1531 else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
1535 printf("The Underflow threshold is %.17e, %s\n", UfThold,
1537 printf("calculation may suffer larger Relative error than ");
1538 printf("merely roundoff.\n");
1542 if (Y2 <= UfThold) {
1544 BadCond(Defect, "");
1548 BadCond(Serious, "");
1551 printf("Range is too narrow; U1^%d Underflows.\n", I);
1553 /*=============================================*/
1556 #include "paranoia.h"
1560 /*=============================================*/
1561 Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
1563 printf("Since underflow occurs below the threshold\n");
1564 printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
1565 printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n",
1567 printf("actually calculating yields:");
1568 if (setjmp(ovfl_buf)) {
1570 BadCond(Serious, "trap on underflow.\n");
1574 V9 = POW(HInvrse, Y2);
1576 printf(" %.17e .\n", V9);
1577 if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
1578 BadCond(Serious, "this is not between 0 and underflow\n");
1579 printf(" threshold = %.17e .\n", UfThold);
1581 else if (! (V9 > UfThold * (One + E9)))
1582 printf("This computed value is O.K.\n");
1584 BadCond(Defect, "this is not between 0 and underflow\n");
1585 printf(" threshold = %.17e .\n", UfThold);
1588 /*=============================================*/
1590 /*=============================================*/
1592 /* ...calculate Exp2 == exp(2) == 7.389056099... */
1606 Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
1611 printf("X: %.30f BInverse: %.30f Radix: %.30f Exp2: %.30f F9: %.30f U1 %.30f U2 %.30f F9: %.30f\n", X, BInvrse, Radix, Exp2, F9, U1, U2, F9);
1612 printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
1616 Z = (X + One) / (Z - (One - BInvrse));
1617 Q = POW(X, Z) - Exp2;
1618 if (FABS(Q) > TwoForty * U2) {
1620 V9 = (X - BInvrse) - (One - BInvrse);
1621 BadCond(Defect, "Calculated");
1622 printf(" %.17e for\n", POW(X,Z));
1623 printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
1624 printf("\tdiffers from correct value by %.17e .\n", Q);
1625 printf("\tThis much error may spoil financial\n");
1626 printf("\tcalculations involving tiny interest rates.\n");
1630 Z = (Y - X) * Two + Y;
1633 Z = One + (X - F9)*(X - F9);
1634 if (Z > One && I < NoTrials) I++;
1638 printf("Accuracy seems adequate.\n");
1650 /*=============================================*/
1652 /*=============================================*/
1653 printf("Testing powers Z^Q at four nearly extreme values.\n");
1656 Q = FLOOR(Half - LOG(C) / LOG(A1));
1666 if (Z < One) Break = True;
1668 } while ( ! (Break));
1670 if (N == 0) printf(" ... no discrepancies found.\n");
1673 /*=============================================*/
1675 /*=============================================*/
1677 printf("Searching for Overflow threshold:\n");
1678 printf("This may generate an error.\n");
1682 if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
1692 printf("Can `Z = -Y' overflow?\n");
1693 printf("Trying it on Y = %.17e .\n", Y);
1696 if (V - Y == V + V0) printf("Seems O.K.\n");
1699 BadCond(Flaw, "-(-Y) differs from Y.\n");
1702 BadCond(Serious, "");
1703 printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
1706 Y = V * (HInvrse * U2 - HInvrse);
1707 Z = Y + ((One - HInvrse) * U2) * V;
1710 if (V0 - V < V0) V = V0;
1713 V = Y * (HInvrse * U2 - HInvrse);
1714 V = V + ((One - HInvrse) * U2) * Y;
1716 printf("Overflow threshold is V = %.17e .\n", V);
1717 if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
1718 else printf("There is no saturation value because \
1719 the system traps on overflow.\n");
1721 printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
1723 printf(" nor for V / 1 = %.17e .\n", V9);
1724 printf("Any overflow signal separating this * from the one\n");
1725 printf("above is a DEFECT.\n");
1726 /*=============================================*/
1728 /*=============================================*/
1729 if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
1730 BadCond(Failure, "Comparisons involving ");
1731 printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
1734 /*=============================================*/
1736 /*=============================================*/
1738 for(Indx = 1; Indx <= 3; ++Indx) {
1740 case 1: Z = UfThold; break;
1741 case 2: Z = E0; break;
1742 case 3: Z = PseudoZero; break;
1747 if (Y / (One - Radix * E9) < Z
1748 || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
1749 if (V9 > U1) BadCond(Serious, "");
1750 else BadCond(Defect, "");
1751 printf("Comparison alleges that what prints as Z = %.17e\n", Z);
1752 printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
1756 /*=============================================*/
1758 /*=============================================*/
1759 for(Indx = 1; Indx <= 2; ++Indx) {
1760 if (Indx == 1) Z = V;
1763 X = (One - Radix * E9) * V9;
1765 if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
1767 if (X < W) BadCond(Serious, "");
1768 else BadCond(Defect, "");
1769 printf("Comparison alleges that Z = %17e\n", Z);
1770 printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
1773 /*=============================================*/
1776 #include "paranoia.h"
1780 /*=============================================*/
1784 if (X*Y < One || X > Y) {
1785 if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
1786 else BadCond(Flaw, "");
1788 printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
1789 X, "is too far from 1.\n");
1791 /*=============================================*/
1793 /*=============================================*/
1794 for (Indx = 1; Indx <= 5; ++Indx) {
1797 case 2: X = One + U2; break;
1798 case 3: X = V; break;
1799 case 4: X = UfThold; break;
1804 if (setjmp(ovfl_buf))
1805 printf(" X / X traps when X = %g\n", X);
1807 V9 = (Y / X - Half) - Half;
1808 if (V9 == Zero) continue;
1809 if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
1810 else BadCond(Serious, "");
1811 printf(" X / X differs from 1 when X = %.17e\n", X);
1812 printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
1816 /*=============================================*/
1818 /*=============================================*/
1821 printf("What message and/or values does Division by Zero produce?\n") ;
1823 printf("This can interupt your program. You can ");
1824 printf("skip this part if you wish.\n");
1825 printf("Do you wish to compute 1 / 0? ");
1827 read (KEYBOARD, ch, 8);
1828 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1831 printf(" Trying to compute 1 / 0 produces ...");
1832 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
1836 else printf("O.K.\n");
1837 printf("\nDo you wish to compute 0 / 0? ");
1839 read (KEYBOARD, ch, 80);
1840 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1843 printf("\n Trying to compute 0 / 0 produces ...");
1844 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
1848 else printf("O.K.\n");
1850 /*=============================================*/
1852 /*=============================================*/
1856 static char *msg[] = {
1857 "FAILUREs encountered =",
1858 "SERIOUS DEFECTs discovered =",
1859 "DEFECTs discovered =",
1860 "FLAWs discovered =" };
1862 for(i = 0; i < 4; i++) if (ErrCnt[i])
1863 printf("The number of %-29s %d.\n",
1867 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
1868 + ErrCnt[Flaw]) > 0) {
1869 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
1870 Defect] == 0) && (ErrCnt[Flaw] > 0)) {
1871 printf("The arithmetic diagnosed seems ");
1872 printf("Satisfactory though flawed.\n");
1874 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
1875 && ( ErrCnt[Defect] > 0)) {
1876 printf("The arithmetic diagnosed may be Acceptable\n");
1877 printf("despite inconvenient Defects.\n");
1879 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
1880 printf("The arithmetic diagnosed has ");
1881 printf("unacceptable Serious Defects.\n");
1883 if (ErrCnt[Failure] > 0) {
1884 printf("Potentially fatal FAILURE may have spoiled this");
1885 printf(" program's subsequent diagnoses.\n");
1889 printf("No failures, defects nor flaws have been discovered.\n");
1890 if (! ((RMult == Rounded) && (RDiv == Rounded)
1891 && (RAddSub == Rounded) && (RSqrt == Rounded)))
1892 printf("The arithmetic diagnosed seems Satisfactory.\n");
1894 if (StickyBit >= One &&
1895 (Radix - Two) * (Radix - Nine - One) == Zero) {
1896 printf("Rounding appears to conform to ");
1897 printf("the proposed IEEE standard P");
1898 if ((Radix == Two) &&
1899 ((Precision - Four * Three * Two) *
1900 ( Precision - TwentySeven -
1901 TwentySeven + One) == Zero))
1904 if (IEEE) printf(".\n");
1906 printf(",\nexcept for possibly Double Rounding");
1907 printf(" during Gradual Underflow.\n");
1910 printf("The arithmetic diagnosed appears to be Excellent!\n");
1914 printf("\nA total of %d floating point exceptions were registered.\n",
1916 printf("END OF TEST.\n");
1921 #include "paranoia.h"
1929 { return X >= (FLOAT)0. ? (FLOAT)1.0 : (FLOAT)-1.0; }
1937 printf("\nTo continue, press RETURN");
1939 read(KEYBOARD, ch, 8);
1941 printf("\nDiagnosis resumes after milestone Number %d", Milestone);
1942 printf(" Page: %d\n\n", PageNo);
1948 TstCond (INT K, INT Valid, CHARP T)
1953 { if (! Valid) { BadCond(K,T); printf(".\n"); } }
1956 BadCond(INT K, CHARP T)
1962 static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
1964 ErrCnt [K] = ErrCnt [K] + 1;
1965 printf("%s: %s", msg[K], T);
1972 X = (Random1 + Random9)^5
1973 Random1 = X - FLOOR(X) + 0.000005 * X;
1974 and returns the new value of Random1
1979 X = Random1 + Random9;
1984 Random1 = Y + X * (FLOAT)0.000005;
1989 SqXMinX (INT ErrKind)
1998 SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
2000 if (SqEr < MinSqEr) MinSqEr = SqEr;
2001 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
2003 BadCond(ErrKind, "\n");
2004 printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
2005 printf("\tinstead of correct value 0 .\n");
2013 X = FLOOR(Half - X / Radix) * Radix + X;
2014 Q = (Q - X * Z) / Radix + X * X * (D / Radix);
2015 Z = Z - Two * X * D;
2026 if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
2029 Y2 = (X2 - Z2) - (Y - Z2);
2030 X2 = X8 / (Y - Half);
2031 X2 = X2 - Half * X2 * X2;
2032 SqEr = (Y2 + Half) + (Half - X2);
2033 if (SqEr < MinSqEr) MinSqEr = SqEr;
2035 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
2044 if (Z == Zero && Q <= Zero)
2045 printf("WARNING: computing\n");
2046 else BadCond(Defect, "computing\n");
2047 printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
2048 printf("\tyielded %.17e;\n", Y);
2049 printf("\twhich compared unequal to correct %.17e ;\n",
2051 printf("\t\tthey differ by %.17e .\n", Y - X);
2053 N = N + 1; /* ... count discrepancies. */
2070 PrintIfNPositive(VOID)
2072 if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
2080 printf("Since comparison denies Z = 0, evaluating ");
2081 printf("(Z + Z) / Z should be safe.\n");
2083 if (setjmp(ovfl_buf)) goto very_serious;
2085 printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
2087 if (FABS(Q9 - Two) < Radix * U2) {
2088 printf("This is O.K., provided Over/Underflow");
2089 printf(" has NOT just been signaled.\n");
2092 if ((Q9 < One) || (Q9 > Two)) {
2095 ErrCnt [Serious] = ErrCnt [Serious] + 1;
2096 printf("This is a VERY SERIOUS DEFECT!\n");
2100 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2101 printf("This is a DEFECT!\n");
2110 if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
2115 BadCond(Defect, "What prints as Z = ");
2116 printf("%.17e\n\tcompares different from ", Z);
2117 if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
2118 if (! ((Z == Random2)
2119 || (Random2 == Random1)))
2120 printf("1 * Z == %g\n", Random2);
2121 if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
2122 if (Random2 != Random1) {
2123 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2124 BadCond(Defect, "Multiplication does not commute!\n");
2125 printf("\tComparison alleges that 1 * Z = %.17e\n",
2127 printf("\tdiffers from Z * 1 = %.17e\n", Random1);
2140 printf("%s test appears to be inconsistent...\n", s);
2141 printf(" PLEASE NOTIFY KARPINKSI!\n");
2145 #include "paranoia.h"
2153 { while(*s) printf("%s\n", *s++); }
2158 static char *instr[] = {
2159 "Lest this program stop prematurely, i.e. before displaying\n",
2160 " `END OF TEST',\n",
2161 "try to persuade the computer NOT to terminate execution when an",
2162 "error like Over/Underflow or Division by Zero occurs, but rather",
2163 "to persevere with a surrogate value after, perhaps, displaying some",
2164 "warning. If persuasion avails naught, don't despair but run this",
2165 "program anyway to see how many milestones it passes, and then",
2166 "amend it to make further progress.\n",
2167 "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
2176 static char *head[] = {
2177 "Users are invited to help debug and augment this program so it will",
2178 "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
2179 "Please send suggestions and interesting results to",
2180 "\tRichard Karpinski",
2181 "\tComputer Center U-76",
2182 "\tUniversity of California",
2183 "\tSan Francisco, CA 94143-0704, USA\n",
2184 "In doing so, please include the following information:",
2186 "\tPrecision:\tsingle;",
2188 "\tPrecision:\tdouble;",
2190 "\tVersion:\t10 February 1989;",
2193 "\tOptimization level:\n",
2194 "\tOther relevant compiler options:",
2201 Characteristics(VOID)
2203 static char *chars[] = {
2204 "Running this program should reveal these characteristics:",
2205 " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
2206 " Precision = number of significant digits carried.",
2207 " U2 = Radix/Radix^Precision = One Ulp",
2208 "\t(OneUlpnit in the Last Place) of 1.000xxx .",
2209 " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
2210 " Adequacy of guard digits for Mult., Div. and Subt.",
2211 " Whether arithmetic is chopped, correctly rounded, or something else",
2212 "\tfor Mult., Div., Add/Subt. and Sqrt.",
2213 " Whether a Sticky Bit used correctly for rounding.",
2214 " UnderflowThreshold = an underflow threshold.",
2215 " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
2216 " V = an overflow threshold, roughly.",
2217 " V0 tells, roughly, whether Infinity is represented.",
2218 " Comparisions are checked for consistency with subtraction",
2219 "\tand for contamination with pseudo-zeros.",
2220 " Sqrt is tested. Y^X is not tested.",
2221 " Extra-precise subexpressions are revealed but NOT YET tested.",
2222 " Decimal-Binary conversion is NOT YET tested for accuracy.",
2231 /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
2232 with further massaging by David M. Gay. */
2234 static char *hist[] = {
2235 "The program attempts to discriminate among",
2236 " FLAWs, like lack of a sticky bit,",
2237 " Serious DEFECTs, like lack of a guard digit, and",
2238 " FAILUREs, like 2+2 == 5 .",
2239 "Failures may confound subsequent diagnoses.\n",
2240 "The diagnostic capabilities of this program go beyond an earlier",
2241 "program called `MACHAR', which can be found at the end of the",
2242 "book `Software Manual for the Elementary Functions' (1980) by",
2243 "W. J. Cody and W. Waite. Although both programs try to discover",
2244 "the Radix, Precision and range (over/underflow thresholds)",
2245 "of the arithmetic, this program tries to cope with a wider variety",
2246 "of pathologies, and to say how well the arithmetic is implemented.",
2247 "\nThe program is based upon a conventional radix representation for",
2248 "floating-point numbers, but also allows logarithmic encoding",
2249 "as used by certain early WANG machines.\n",
2250 "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
2251 "see source comments for more history.",