2 !!DESCRIPTION!! a wellknown floating point test
3 !!ORIGIN!! LCC 4.1 Testsuite
4 !!LICENCE!! own, freely distributeable for non-profit. read CPYRIGHT.LCC
13 printf("NO_FLOATS\n\r");
21 /* A C version of Kahan's Floating Point Test "Paranoia"
23 Thos Sumner, UCSF, Feb. 1985
24 David Gay, BTL, Jan. 1986
26 This is a rewrite from the Pascal version by
28 B. A. Wichmann, 18 Jan. 1985
30 (and does NOT exhibit good C programming style).
32 (C) Apr 19 1983 in BASIC version by:
33 Professor W. M. Kahan,
35 Electrical Engineering & Computer Science Dept.
36 University of California
37 Berkeley, California 94720
40 converted to Pascal by:
42 National Physical Laboratory
49 David M. Gay and Thos Sumner
50 AT&T Bell Labs Computer Center, Rm. U-76
51 600 Mountain Avenue University of California
52 Murray Hill, NJ 07974 San Francisco, CA 94143
55 with simultaneous corrections to the Pascal source (reflected
56 in the Pascal source available over netlib).
57 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
59 Reports of results on various systems from all the versions
60 of Paranoia are being collected by Richard Karpinski at the
61 same address as Thos Sumner. This includes sample outputs,
62 bug reports, and criticisms.
64 You may copy this program freely if you acknowledge its source.
65 Comments on the Pascal version to NPL, please.
67 The C version catches signals from floating-point exceptions.
68 If signal(SIGFPE,...) is unavailable in your environment, you may
69 #define NOSIGNAL to comment out the invocations of signal.
71 This source file is too big for some C compilers, but may be split
72 into pieces. Comments containing "SPLIT" suggest convenient places
73 for this splitting. At the end of these comments is an "ed script"
74 (for the UNIX(tm) editor ed) that will do this splitting.
76 By #defining Single when you compile this source, you may obtain
77 a single-precision C version of Paranoia.
79 The following is from the introductory commentary from Wichmann's work:
81 The BASIC program of Kahan is written in Microsoft BASIC using many
82 facilities which have no exact analogy in Pascal. The Pascal
83 version below cannot therefore be exactly the same. Rather than be
84 a minimal transcription of the BASIC program, the Pascal coding
85 follows the conventional style of block-structured languages. Hence
86 the Pascal version could be useful in producing versions in other
89 Rather than use identifiers of minimal length (which therefore have
90 little mnemonic significance), the Pascal version uses meaningful
91 identifiers as follows [Note: A few changes have been made for C]:
93 BASIC C BASIC C BASIC C
96 A1 AInverse J0 NoErrors T
97 B Radix [Failure] T0 Underflow
98 B1 BInverse J1 NoErrors T2 ThirtyTwo
99 B2 RadixD2 [SeriousDefect] T5 OneAndHalf
100 B9 BMinusU2 J2 NoErrors T7 TwentySeven
101 C [Defect] T8 TwoForty
102 C1 CInverse J3 NoErrors U OneUlp
103 D [Flaw] U0 UnderflowThreshold
114 F1 MinusOne O5 Five X9 Random1
120 G2 GDiv Q9 Z0 PseudoZero
123 H1 HInverse R2 RDiv Z9
130 All the variables in BASIC are true variables and in consequence,
131 the program is more difficult to follow since the "constants" must
132 be determined (the glossary is very helpful). The Pascal version
133 uses Real constants, but checks are added to ensure that the values
134 are correctly converted by the compiler.
136 The major textual change to the Pascal version apart from the
137 identifiersis that named procedures are used, inserting parameters
138 wherehelpful. New procedures are also introduced. The
139 correspondence is as follows:
145 170- 250 Instructions
147 480- 670 Characteristics
151 4040-4080 DoesYequalX
152 4090-4110 PrintIfNPositive
153 4640-4850 TestPartialUnderflow
155 =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
157 Below is an "ed script" that splits para.c into 10 files
158 of the form part[1-8].c, subs.c, and msgs.c, plus a header
159 file, paranoia.h, that these files require.
216 /Computed constants/,$d
218 1,$s/^FLOAT/extern &/
221 /^Guard/,/^Round/s/^/extern /
222 /^jmp_buf/s/^/extern /
223 /^Sig_type/s/^/extern /
225 extern void sigfpe();/
237 extern double fabs(), floor(), log(), pow(), sqrt();
241 #define FABS(x) (float)fabs((double)(x))
242 #define FLOOR(x) (float)floor((double)(x))
243 #define LOG(x) (float)log((double)(x))
244 #define POW(x,y) (float)pow((double)(x),(double)(y))
245 #define SQRT(x) (float)sqrt((double)(x))
248 #define FABS(x) fabs(x)
249 #define FLOOR(x) floor(x)
250 #define LOG(x) log(x)
251 #define POW(x,y) pow(x,y)
252 #define SQRT(x) sqrt(x)
256 typedef void (*Sig_type)();
261 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
262 FLOAT Sign(), Random();
264 /*Small floating point constants.*/
274 FLOAT TwentySeven = 27.0;
275 FLOAT ThirtyTwo = 32.0;
276 FLOAT TwoForty = 240.0;
277 FLOAT MinusOne = -1.0;
278 FLOAT OneAndHalf = 1.5;
279 /*Integer constants*/
280 int NoTrials = 20; /*Number of tests for commutativity. */
284 /* Definitions for declared types
286 Rounding == (Chopped, Rounded, Other);
287 Message == packed array [1..40] of char;
288 Class == (Flaw, Defect, Serious, Failure);
299 typedef int Guard, Rounding, Class;
300 typedef char Message;
302 /* Declarations of Variables */
308 FLOAT E0, E1, Exp2, E3, MinSqEr;
309 FLOAT SqEr, MaxSqEr, E9;
319 FLOAT T, Underflow, S;
320 FLOAT OneUlp, UfThold, U1, U2;
323 FLOAT X, X1, X2, X8, Random1;
324 FLOAT Y, Y1, Y2, Random2;
325 FLOAT Z, PseudoZero, Z1, Z2, Z9;
331 Guard GMult, GDiv, GAddSub;
332 Rounding RMult, RDiv, RAddSub, RSqrt;
333 int Break, Done, NotMonot, Monot, Anomaly, IEEE,
335 /* Computed constants. */
336 /*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
337 /*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
339 /* floating point exception receiver */
344 printf("\n* * * FLOATING-POINT ERROR * * *\n");
348 signal(SIGFPE, sigsave);
351 longjmp(ovfl_buf, 1);
360 ieee_flags("set", "precision", "double", &out);
362 /* First two assignments use integer right-hand sides. */
370 Nine = Three * Three;
371 TwentySeven = Nine * Three;
372 ThirtyTwo = Four * Eight;
373 TwoForty = Four * Five * Three * Four;
376 OneAndHalf = One + Half;
382 /*=============================================*/
384 /*=============================================*/
386 signal(SIGFPE, sigfpe);
396 /*=============================================*/
398 /*=============================================*/
399 printf("Program is now RUNNING tests on small integers:\n");
401 TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
402 && (One > Zero) && (One + One == Two),
403 "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
406 ErrCnt[Failure] = ErrCnt[Failure] + 1;
407 printf("Comparison alleges that -0.0 is Non-zero!\n");
412 TstCond (Failure, (Three == Two + One) && (Four == Three + One)
413 && (Four + Two * (- Two) == Zero)
414 && (Four - Three - One == Zero),
415 "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
416 TstCond (Failure, (MinusOne == (0 - One))
417 && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
418 && (MinusOne + FABS(One) == Zero)
419 && (MinusOne + MinusOne * MinusOne == Zero),
420 "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
421 TstCond (Failure, Half + MinusOne + Half == Zero,
422 "1/2 + (-1) + 1/2 != 0");
423 /*=============================================*/
433 #include "paranoia.h"
437 /*=============================================*/
438 TstCond (Failure, (Nine == Three * Three)
439 && (TwentySeven == Nine * Three) && (Eight == Four + Four)
440 && (ThirtyTwo == Eight * Four)
441 && (ThirtyTwo - TwentySeven - Four - One == Zero),
442 "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
443 TstCond (Failure, (Five == Four + One) &&
444 (TwoForty == Four * Five * Three * Four)
445 && (TwoForty / Three - Four * Four * Five == Zero)
446 && ( TwoForty / Four - Five * Three * Four == Zero)
447 && ( TwoForty / Five - Four * Three * Four == Zero),
448 "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
449 if (ErrCnt[Failure] == 0) {
450 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
453 printf("Searching for Radix and Precision.\n");
460 } while (MinusOne + FABS(Y) < Zero);
461 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
468 } while ( Radix == Zero);
469 if (Radix < Two) Radix = One;
470 printf("Radix = %f .\n", Radix);
474 Precision = Precision + One;
477 } while ((Y - W) == One);
479 /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
483 printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
484 printf("Recalculating radix and precision\n ");
499 /*... now X = (unknown no.) ulps of 1+...*/
502 Y = Half * U2 + ThirtyTwo * U2 * U2;
505 } while ( ! ((U2 <= X) || (X <= Zero)));
507 /*... now U2 == 1 ulp of 1 + ... */
515 /*... now X == (unknown no.) ulps of 1 -... */
518 Y = Half * U1 + ThirtyTwo * U1 * U1;
523 } while ( ! ((U1 <= X) || (X <= Zero)));
524 /*... now U1 == 1 ulp of 1 - ... */
525 if (U1 == E1) printf("confirms closest relative separation U1 .\n");
526 else printf("gets better closest relative separation U1 = %.7e .\n", U1);
528 F9 = (Half - U1) + Half;
529 Radix = FLOOR(0.01 + U2 / U1);
530 if (Radix == E0) printf("Radix confirmed.\n");
531 else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
532 TstCond (Defect, Radix <= Eight + Eight,
533 "Radix is too big: roundoff problems");
534 TstCond (Flaw, (Radix == Two) || (Radix == 10)
535 || (Radix == One), "Radix is not as good as 2 or 10");
536 /*=============================================*/
538 /*=============================================*/
539 TstCond (Failure, F9 - Half < Half,
540 "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
545 TstCond (Failure, (X != One)
546 || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
549 /*=============================================*/
551 /*=============================================*/
552 /*... BMinusU2 = nextafter(Radix, 0) */
553 BMinusU2 = Radix - One;
554 BMinusU2 = (BMinusU2 - U2) + One;
555 /* Purify Integers */
557 X = - TwoForty * LOG(U1) / LOG(Radix);
559 if (FABS(X - Y) * Four < One) X = Y;
560 Precision = X / TwoForty;
561 Y = FLOOR(Half + Precision);
562 if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
564 if ((Precision != FLOOR(Precision)) || (Radix == One)) {
565 printf("Precision cannot be characterized by an Integer number\n");
566 printf("of significant digits but, by itself, this is a minor flaw.\n");
569 printf("logarithmic encoding has precision characterized solely by U1.\n");
570 else printf("The number of significant digits of the Radix is %f .\n",
572 TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
573 "Precision worse than 5 decimal figures ");
574 /*=============================================*/
576 /*=============================================*/
577 /* Test for extra-precise subepressions */
578 X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
581 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
582 } while ( ! ((Z2 <= X) || (X <= Zero)));
583 X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
586 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
587 + One / Two)) + One / Two;
588 } while ( ! ((Z1 <= Z) || (Z <= Zero)));
592 Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
594 } while ( ! ((Y1 <= Y) || (Y <= Zero)));
596 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
597 } while ( ! ((X1 <= X) || (X <= Zero)));
598 if ((X1 != Y1) || (X1 != Z1)) {
599 BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
600 printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
601 printf("are symptoms of inconsistencies introduced\n");
602 printf("by extra-precise evaluation of arithmetic subexpressions.\n");
603 notify("Possibly some part of this");
604 if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
605 "That feature is not tested further by this program.\n") ;
608 if ((Z1 != U1) || (Z2 != U2)) {
609 if ((Z1 >= U1) || (Z2 >= U2)) {
610 BadCond(Failure, "");
612 printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
613 printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
616 if ((Z1 <= Zero) || (Z2 <= Zero)) {
617 printf("Because of unusual Radix = %f", Radix);
618 printf(", or exact rational arithmetic a result\n");
619 printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
620 notify("of an\nextra-precision");
622 if (Z1 != Z2 || Z1 > Zero) {
627 printf("Some subexpressions appear to be calculated extra\n");
628 printf("precisely with about %g extra B-digits, i.e.\n",
630 printf("roughly %g extra significant decimals.\n",
633 printf("That feature is not tested further by this program.\n");
638 /*=============================================*/
641 #include "paranoia.h"
645 /*=============================================*/
647 X = W / (Radix * Radix);
652 TstCond (Failure, X == U2,
653 "Subtraction is not normalized X=Y,X+Z != Y+Z!");
655 "Subtraction appears to be normalized, as it should be.");
657 printf("\nChecking for guard digit in *, /, and -.\n");
670 X = X * (Radix - One);
671 T = T * (Radix - One);
672 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
675 TstCond (Serious, False,
676 "* lacks a Guard Digit, so 1*X != X");
680 Y = FABS((X + Z) - X * X) - U2;
682 Z = FABS((X - U2) - X * X) - U1;
683 TstCond (Failure, (Y <= Zero)
684 && (Z <= Zero), "* gets too many final digits wrong.\n");
692 T = Nine / TwentySeven;
694 TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
695 "Division lacks a Guard Digit, so error can exceed 1 ulp\nor 1/3 and 3/9 and 9/27 may disagree");
702 if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
705 TstCond (Serious, False,
706 "Division lacks a Guard Digit, so X/1 != X");
708 X = One / (One + U2);
710 TstCond (Serious, Y < Zero,
711 "Computed value of 1/1.000..1 >= 1");
713 Y = One + Radix * U2;
717 StickyBit = T / Radix;
720 TstCond (Failure, X == Zero && Y == Zero,
721 "* and/or / gets too many last digits wrong");
726 Z = Radix - BMinusU2;
728 if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
731 TstCond (Serious, False,
732 "- lacks Guard Digit, so cancellation is obscured");
734 if (F9 != One && F9 - One >= Zero) {
735 BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
736 printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
737 printf(" such precautions against division by zero as\n");
738 printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
740 if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
741 " *, /, and - appear to have guard digits, as they should.\n");
742 /*=============================================*/
744 /*=============================================*/
746 printf("Checking rounding on multiply, divide and add/subtract.\n");
750 RadixD2 = Radix / Two;
757 AInvrse = AInvrse / A1;
758 } while ( ! (FLOOR(AInvrse) != AInvrse));
759 Done = (X == One) || (A1 > Three);
760 if (! Done) A1 = Nine + One;
762 if (X == One) A1 = Radix;
769 TstCond (Failure, Z == Half,
770 "X * (1/X) differs from 1");
787 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
788 X = (OneAndHalf + U2) * Y2;
789 Y = OneAndHalf - U2 - U2;
790 Z = OneAndHalf + U2 + U2;
791 T = (OneAndHalf - U2) * Y1;
796 Y = (U2 - Y) + StickyBit;
797 Z = S - (Z + U2 + U2);
798 StickyBit = (Y2 + U2) * Y1;
800 StickyBit = StickyBit - Y2;
802 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
803 && ( StickyBit == Zero) && (Y1 == Half)) {
805 printf("Multiplication appears to round correctly.\n");
807 else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
808 && (T < Zero) && (StickyBit + U2 == Zero)
811 printf("Multiplication appears to chop.\n");
813 else printf("* is neither chopped nor correctly rounded.\n");
814 if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
816 else printf("* is neither chopped nor correctly rounded.\n");
817 /*=============================================*/
819 /*=============================================*/
822 Z = OneAndHalf + U2 + U2;
824 T = OneAndHalf - U2 - U2;
830 Z = Z - (OneAndHalf + U2);
831 T = (U2 - OneAndHalf) + T;
832 if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
846 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
847 && (Y2 == Zero) && (Y2 == Zero)
848 && (Y1 - Half == F9 - Half )) {
850 printf("Division appears to round correctly.\n");
851 if (GDiv == No) notify("Division");
853 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
854 && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
856 printf("Division appears to chop.\n");
859 if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
860 BInvrse = One / Radix;
861 TstCond (Failure, (BInvrse * Radix - Half == Half),
862 "Radix * ( 1 / Radix ) differs from 1");
863 /*=============================================*/
866 #include "paranoia.h"
870 /*=============================================*/
871 TstCond (Failure, ((F9 + U1) - Half == Half)
872 && ((BMinusU2 + U2 ) - One == Radix - One),
873 "Incomplete carry-propagation in Addition");
875 Y = One + U2 * (One - U2);
879 if ((X == Zero) && (Y == Zero)) {
881 printf("Add/Subtract appears to be chopped.\n");
883 if (GAddSub == Yes) {
884 X = (Half + U2) * U2;
885 Y = (Half - U2) * U2;
890 if ((X == Zero) && (Y == Zero)) {
891 X = (Half + U2) * U1;
892 Y = (Half - U2) * U1;
897 if ((X == Zero) && (Y == Zero)) {
899 printf("Addition/Subtraction appears to round correctly.\n");
900 if (GAddSub == No) notify("Add/Subtract");
902 else printf("Addition/Subtraction neither rounds nor chops.\n");
904 else printf("Addition/Subtraction neither rounds nor chops.\n");
906 else printf("Addition/Subtraction neither rounds nor chops.\n");
908 X = One + Half * (One + Half);
909 Y = (One + U2) * Half;
913 if (StickyBit != Zero) {
915 BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
918 if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
919 && (RMult == Rounded) && (RDiv == Rounded)
920 && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
921 printf("Checking for sticky bit.\n");
922 X = (Half + U1) * U2;
926 if ((Z - One <= Zero) && (T - One >= U2)) {
929 if ((Z - T >= U2) && (Y - T == Zero)) {
930 X = (Half + U1) * U1;
934 if ((Z - One == Zero) && (T - F9 == Zero)) {
935 Z = (Half - U1) * U1;
938 if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
939 Z = (One + U2) * OneAndHalf;
940 T = (OneAndHalf + U2) - Z + U2;
941 X = One + Half / Radix;
942 Y = One + Radix * U2;
944 if (T == Zero && X + Radix * U2 - Z == Zero) {
948 if ((Y - One == Zero)) StickyBit = S;
957 if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
958 else printf("Sticky bit used incorrectly or not at all.\n");
959 TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
960 RMult == Other || RDiv == Other || RAddSub == Other),
961 "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");
962 /*=============================================*/
964 /*=============================================*/
966 printf("Does Multiplication commute? ");
967 printf("Testing on %d random pairs.\n", NoTrials);
978 } while ( ! ((I > NoTrials) || (Z9 != Zero)));
980 Random1 = One + Half / Three;
981 Random2 = (U2 + U1) + One;
982 Z = Random1 * Random2;
983 Y = Random2 * Random1;
984 Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
985 Three) * ((U2 + U1) + One);
987 if (! ((I == NoTrials) || (Z9 == Zero)))
988 BadCond(Defect, "X * Y == Y * X trial fails.\n");
989 else printf(" No failures found in %d integer pairs.\n", NoTrials);
990 /*=============================================*/
992 /*=============================================*/
993 printf("\nRunning test of square root(x).\n");
994 TstCond (Failure, (Zero == SQRT(Zero))
995 && (- Zero == SQRT(- Zero))
996 && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
1004 OneUlp = BInvrse * U1;
1009 if (J != Zero) Pause();
1010 printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
1014 if ((Radix != One)) do {
1017 } while ( ! ((Y - X >= NoTrials)));
1020 while (I <= NoTrials) {
1023 if (J > Zero) break;
1026 printf("Test for sqrt monotonicity.\n");
1030 Z = Radix + Radix * U2;
1033 while ( ! (NotMonot || Monot)) {
1038 if ((X > Q) || (Q > Z)) NotMonot = True;
1040 Q = FLOOR(Q + Half);
1041 if ((I > 0) || (Radix == Q * Q)) Monot = True;
1043 if (I > 1) Monot = True;
1057 if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
1059 BadCond(Defect, "");
1060 printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
1062 /*=============================================*/
1065 #include "paranoia.h"
1069 /*=============================================*/
1070 MinSqEr = MinSqEr + Half;
1071 MaxSqEr = MaxSqEr - Half;
1072 Y = (SQRT(One + U2) - One) / U2;
1073 SqEr = (Y - One) + U2 / Eight;
1074 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1075 SqEr = Y + U2 / Eight;
1076 if (SqEr < MinSqEr) MinSqEr = SqEr;
1077 Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
1078 SqEr = Y + U1 / Eight;
1079 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1080 SqEr = (Y + One) + U1 / Eight;
1081 if (SqEr < MinSqEr) MinSqEr = SqEr;
1084 for( Indx = 1; Indx <= 3; ++Indx) {
1085 Y = SQRT((X + U1 + X) + F9);
1086 Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
1087 Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
1088 SqEr = (Y + Half) + Z;
1089 if (SqEr < MinSqEr) MinSqEr = SqEr;
1090 SqEr = (Y - Half) + Z;
1091 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1092 if (((Indx == 1) || (Indx == 3)))
1093 X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
1099 /*=============================================*/
1101 /*=============================================*/
1104 RSqrt = Other; /* ~dgh */
1106 printf("Testing whether sqrt is rounded or chopped.\n");
1107 D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
1108 /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
1111 if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
1126 X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
1129 } while ( ! (X1 <= Zero));
1138 if (Y2 >= FourD) Y2 = Y2 - FourD;
1139 } while ( ! (Y >= D));
1141 Q = (X8 + Z * Z) / FourD;
1143 if (Q != FLOOR(Q)) Anomaly = True;
1148 X = X - FLOOR(X / Radix) * Radix;
1153 } while ( ! (Break || (Z1 <= Zero)));
1154 if ((Z1 <= Zero) && (! Break)) Anomaly = True;
1156 if (Z1 > RadixD2) Z1 = Z1 - Radix;
1159 } while ( ! (U2 * D >= F9));
1160 if (D * Radix - D != W - D) Anomaly = True;
1164 Y = D + (One + Z) * Half;
1167 Y = D + (One - Z) * Half + D;
1172 if (D - Z2 != W - Z2) Anomaly = True;
1174 Y = (D - Z2) + (Z2 + (One - Z) * Half);
1175 X = (D - Z2) + (Z2 - Z + Q);
1177 Y = (One + Z) * Half;
1180 if (I == 0) Anomaly = True;
1186 if ((I == 0) || Anomaly) {
1187 BadCond(Failure, "Anomalous arithmetic with Integer < ");
1188 printf("Radix^Precision = %.7e\n", W);
1189 printf(" fails test whether sqrt rounds or chops.\n");
1194 if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
1196 printf("Square root appears to be correctly rounded.\n");
1199 if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
1200 || (MinSqEr + Radix < Half)) SqRWrng = True;
1203 printf("Square root appears to be chopped.\n");
1208 printf("Square root is neither chopped nor correctly rounded.\n");
1209 printf("Observed errors run from %.7e ", MinSqEr - Half);
1210 printf("to %.7e ulps.\n", Half + MaxSqEr);
1211 TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
1212 "sqrt gets too many last digits wrong");
1214 /*=============================================*/
1216 /*=============================================*/
1218 printf("Testing powers Z^i for small Integers Z and i.\n");
1220 /* ... test powers of zero. */
1232 if (Z == MinusOne) Break = True;
1237 /* .. if(-1)^N is invalid, replace MinusOne by One. */
1245 M = FLOOR(Two * LOG(W) / LOG(A1));
1251 if (Z == AInvrse) Break = True;
1253 } while ( ! (Break));
1254 /*=============================================*/
1256 /*=============================================*/
1257 /* Powers of Radix have been tested, */
1258 /* next try a few primes */
1267 } while ( Three * FLOOR(Z / Three) == Z );
1268 } while ( Z < Eight * Three );
1270 printf("Errors like this may invalidate financial calculations\n");
1271 printf("\tinvolving interest rates.\n");
1275 if (N == 0) printf("... no discrepancis found.\n");
1278 /*=============================================*/
1281 #include "paranoia.h"
1285 /*=============================================*/
1286 printf("Seeking Underflow thresholds UfThold and E0.\n");
1288 if (Precision != FLOOR(Precision)) {
1294 } while ( X > Zero);
1298 /* ... D is power of 1/Radix < 1. */
1303 } while ((Y > Z) && (Z + Z > Z));
1310 } while ((Y > Z) && (Z + Z > Z));
1311 if (Radix < Two) HInvrse = Two;
1312 else HInvrse = Radix;
1314 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1318 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1323 } while ((E0 > Z) && (Z + Z > Z));
1335 BadCond(Failure, "multiplication gets too many last digits wrong.\n");
1344 PseudoZero = Underflow * H;
1348 Underflow = PseudoZero;
1349 if (E1 + E1 <= E1) {
1350 Y2 = Underflow * HInvrse;
1353 if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
1355 PseudoZero = PseudoZero * H;
1356 } while ((Underflow > PseudoZero)
1357 && (PseudoZero + PseudoZero > PseudoZero));
1359 /* Comment line 4530 .. 4560 */
1360 if (PseudoZero != Zero) {
1363 /* ... Test PseudoZero for "phoney- zero" violates */
1364 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1366 if (PseudoZero <= Zero) {
1367 BadCond(Failure, "Positive expressions can underflow to an\n");
1368 printf("allegedly negative value\n");
1369 printf("PseudoZero that prints out as: %g .\n", PseudoZero);
1372 printf("But -PseudoZero, which should be\n");
1373 printf("positive, isn't; it prints out as %g .\n", X);
1377 BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
1378 printf("value PseudoZero that prints out as %g .\n", PseudoZero);
1382 /*=============================================*/
1384 /*=============================================*/
1385 if (CInvrse * Y > CInvrse * Y1) {
1389 if (! ((E1 == Zero) || (E1 == E0))) {
1390 BadCond(Defect, "");
1392 printf("Products underflow at a higher");
1393 printf(" threshold than differences.\n");
1394 if (PseudoZero == Zero)
1398 printf("Difference underflows at a higher");
1399 printf(" threshold than products.\n");
1402 printf("Smallest strictly positive number found is E0 = %g .\n", E0);
1406 if (N == 1) Underflow = Y;
1408 if (E1 == Zero) I = 3;
1409 if (UfThold == Zero) I = I - 2;
1413 UfThold = Underflow;
1414 if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
1416 BadCond(Failure, "Either accuracy deteriorates as numbers\n");
1417 printf("approach a threshold = %.17e\n", UfThold);;
1418 printf(" coming down from %.17e\n", C);
1419 printf(" or else multiplication gets too many last digits wrong.\n");
1425 BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
1426 printf("Q == Y while denying that |Q - Y| == 0; these values\n");
1427 printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
1428 printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
1437 if ((Q == UfThold) && (E1 == E0)
1438 && (FABS( UfThold - E1 / E9) <= E1)) {
1440 printf("Underflow is gradual; it incurs Absolute Error =\n");
1441 printf("(roundoff in UfThold) < E0.\n");
1443 Y = Y * (OneAndHalf + U2);
1444 X = CInvrse * (One + U2);
1452 if (setjmp(ovfl_buf)) {
1453 printf("Underflow / UfThold failed!\n");
1456 else R = SQRT(Underflow / UfThold);
1460 X = Z * (One + R * H * (One + H));
1464 X = Z * (One + H * H * (One + H));
1466 if (! ((X == Z) || (X - Z != Zero))) {
1468 printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
1470 printf("yet X - Z yields %.17e .\n", Z9);
1471 printf(" Should this NOT signal Underflow, ");
1472 printf("this is a SERIOUS DEFECT\nthat causes ");
1473 printf("confusion when innocent statements like\n");;
1474 printf(" if (X == Z) ... else");
1475 printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
1476 printf("encounter Division by Zero although actually\n");
1478 if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
1479 else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
1483 printf("The Underflow threshold is %.17e, %s\n", UfThold,
1485 printf("calculation may suffer larger Relative error than ");
1486 printf("merely roundoff.\n");
1490 if (Y2 <= UfThold) {
1492 BadCond(Defect, "");
1496 BadCond(Serious, "");
1499 printf("Range is too narrow; U1^%d Underflows.\n", I);
1501 /*=============================================*/
1504 #include "paranoia.h"
1508 /*=============================================*/
1509 Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
1511 printf("Since underflow occurs below the threshold\n");
1512 printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
1513 printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
1514 V9 = POW(HInvrse, Y2);
1515 printf("actually calculating yields: %.17e .\n", V9);
1516 if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
1517 BadCond(Serious, "this is not between 0 and underflow\n");
1518 printf(" threshold = %.17e .\n", UfThold);
1520 else if (! (V9 > UfThold * (One + E9)))
1521 printf("This computed value is O.K.\n");
1523 BadCond(Defect, "this is not between 0 and underflow\n");
1524 printf(" threshold = %.17e .\n", UfThold);
1526 /*=============================================*/
1528 /*=============================================*/
1530 /* ...calculate Exp2 == exp(2) == 7.389056099... */
1544 Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
1549 printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
1553 Z = (X + One) / (Z - (One - BInvrse));
1554 Q = POW(X, Z) - Exp2;
1555 if (FABS(Q) > TwoForty * U2) {
1557 V9 = (X - BInvrse) - (One - BInvrse);
1558 BadCond(Defect, "Calculated");
1559 printf(" %.17e for\n", POW(X,Z));
1560 printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
1561 printf("\tdiffers from correct value by %.17e .\n", Q);
1562 printf("\tThis much error may spoil financial\n");
1563 printf("\tcalculations involving tiny interest rates.\n");
1567 Z = (Y - X) * Two + Y;
1570 Z = One + (X - F9)*(X - F9);
1571 if (Z > One && I < NoTrials) I++;
1575 printf("Accuracy seems adequate.\n");
1587 /*=============================================*/
1589 /*=============================================*/
1590 printf("Testing powers Z^Q at four nearly extreme values.\n");
1593 Q = FLOOR(Half - LOG(C) / LOG(A1));
1603 if (Z < One) Break = True;
1605 } while ( ! (Break));
1607 if (N == 0) printf(" ... no discrepancies found.\n");
1610 /*=============================================*/
1612 /*=============================================*/
1614 printf("Searching for Overflow threshold:\n");
1615 printf("This may generate an error.\n");
1619 if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
1629 printf("Can `Z = -Y' overflow?\n");
1630 printf("Trying it on Y = %.17e .\n", Y);
1633 if (V - Y == V + V0) printf("Seems O.K.\n");
1636 BadCond(Flaw, "-(-Y) differs from Y.\n");
1639 BadCond(Serious, "");
1640 printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
1643 Y = V * (HInvrse * U2 - HInvrse);
1644 Z = Y + ((One - HInvrse) * U2) * V;
1647 if (V0 - V < V0) V = V0;
1650 V = Y * (HInvrse * U2 - HInvrse);
1651 V = V + ((One - HInvrse) * U2) * Y;
1653 printf("Overflow threshold is V = %.17e .\n", V);
1654 if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
1655 else printf("There is no saturation value because the system traps on overflow.\n");
1657 printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
1659 printf(" nor for V / 1 = %.17e .\n", V9);
1660 printf("Any overflow signal separating this * from the one\n");
1661 printf("above is a DEFECT.\n");
1662 /*=============================================*/
1664 /*=============================================*/
1665 if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
1666 BadCond(Failure, "Comparisons involving ");
1667 printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
1670 /*=============================================*/
1672 /*=============================================*/
1674 for(Indx = 1; Indx <= 3; ++Indx) {
1676 case 1: Z = UfThold; break;
1677 case 2: Z = E0; break;
1678 case 3: Z = PseudoZero; break;
1683 if (Y / (One - Radix * E9) < Z
1684 || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
1685 if (V9 > U1) BadCond(Serious, "");
1686 else BadCond(Defect, "");
1687 printf("Comparison alleges that what prints as Z = %.17e\n", Z);
1688 printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
1692 /*=============================================*/
1694 /*=============================================*/
1695 for(Indx = 1; Indx <= 2; ++Indx) {
1696 if (Indx == 1) Z = V;
1699 X = (One - Radix * E9) * V9;
1701 if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
1703 if (X < W) BadCond(Serious, "");
1704 else BadCond(Defect, "");
1705 printf("Comparison alleges that Z = %17e\n", Z);
1706 printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
1709 /*=============================================*/
1712 #include "paranoia.h"
1716 /*=============================================*/
1720 if (X*Y < One || X > Y) {
1721 if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
1722 else BadCond(Flaw, "");
1724 printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
1725 X, "is too far from 1.\n");
1727 /*=============================================*/
1729 /*=============================================*/
1730 for (Indx = 1; Indx <= 5; ++Indx) {
1733 case 2: X = One + U2; break;
1734 case 3: X = V; break;
1735 case 4: X = UfThold; break;
1740 if (setjmp(ovfl_buf))
1741 printf(" X / X traps when X = %g\n", X);
1743 V9 = (Y / X - Half) - Half;
1744 if (V9 == Zero) continue;
1745 if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
1746 else BadCond(Serious, "");
1747 printf(" X / X differs from 1 when X = %.17e\n", X);
1748 printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
1752 /*=============================================*/
1754 /*=============================================*/
1757 printf("What message and/or values does Division by Zero produce?\n") ;
1759 printf("This can interupt your program. You can ");
1760 printf("skip this part if you wish.\n");
1761 printf("Do you wish to compute 1 / 0? ");
1763 read (KEYBOARD, ch, 8);
1764 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1767 printf(" Trying to compute 1 / 0 produces ...");
1768 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
1772 else printf("O.K.\n");
1773 printf("\nDo you wish to compute 0 / 0? ");
1775 read (KEYBOARD, ch, 80);
1776 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1779 printf("\n Trying to compute 0 / 0 produces ...");
1780 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
1784 else printf("O.K.\n");
1786 /*=============================================*/
1788 /*=============================================*/
1792 static char *msg[] = {
1793 "FAILUREs encountered =",
1794 "SERIOUS DEFECTs discovered =",
1795 "DEFECTs discovered =",
1796 "FLAWs discovered =" };
1798 for(i = 0; i < 4; i++) if (ErrCnt[i])
1799 printf("The number of %-29s %d.\n",
1803 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
1804 + ErrCnt[Flaw]) > 0) {
1805 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
1806 Defect] == 0) && (ErrCnt[Flaw] > 0)) {
1807 printf("The arithmetic diagnosed seems ");
1808 printf("Satisfactory though flawed.\n");
1810 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
1811 && ( ErrCnt[Defect] > 0)) {
1812 printf("The arithmetic diagnosed may be Acceptable\n");
1813 printf("despite inconvenient Defects.\n");
1815 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
1816 printf("The arithmetic diagnosed has ");
1817 printf("unacceptable Serious Defects.\n");
1819 if (ErrCnt[Failure] > 0) {
1820 printf("Potentially fatal FAILURE may have spoiled this");
1821 printf(" program's subsequent diagnoses.\n");
1825 printf("No failures, defects nor flaws have been discovered.\n");
1826 if (! ((RMult == Rounded) && (RDiv == Rounded)
1827 && (RAddSub == Rounded) && (RSqrt == Rounded)))
1828 printf("The arithmetic diagnosed seems Satisfactory.\n");
1830 if (StickyBit >= One &&
1831 (Radix - Two) * (Radix - Nine - One) == Zero) {
1832 printf("Rounding appears to conform to ");
1833 printf("the proposed IEEE standard P");
1834 if ((Radix == Two) &&
1835 ((Precision - Four * Three * Two) *
1836 ( Precision - TwentySeven -
1837 TwentySeven + One) == Zero))
1840 if (IEEE) printf(".\n");
1842 printf(",\nexcept for possibly Double Rounding");
1843 printf(" during Gradual Underflow.\n");
1846 printf("The arithmetic diagnosed appears to be Excellent!\n");
1850 printf("\nA total of %d floating point exceptions were registered.\n",
1852 printf("END OF TEST.\n");
1857 #include "paranoia.h"
1864 { return X >= 0. ? 1.0 : -1.0; }
1873 printf("\nTo continue, press RETURN");
1875 read(KEYBOARD, ch, 8);
1877 printf("\nDiagnosis resumes after milestone Number %d", Milestone);
1878 printf(" Page: %d\n\n", PageNo);
1885 TstCond (K, Valid, T)
1888 { if (! Valid) { BadCond(K,T); printf(".\n"); } }
1894 static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
1896 ErrCnt [K] = ErrCnt [K] + 1;
1897 printf("%s: %s", msg[K], T);
1902 X = (Random1 + Random9)^5
1903 Random1 = X - FLOOR(X) + 0.000005 * X;
1904 and returns the new value of Random1
1911 X = Random1 + Random9;
1916 Random1 = Y + X * 0.000005;
1929 SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
1931 if (SqEr < MinSqEr) MinSqEr = SqEr;
1932 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1934 BadCond(ErrKind, "\n");
1935 printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
1936 printf("\tinstead of correct value 0 .\n");
1945 X = FLOOR(Half - X / Radix) * Radix + X;
1946 Q = (Q - X * Z) / Radix + X * X * (D / Radix);
1947 Z = Z - Two * X * D;
1959 if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
1962 Y2 = (X2 - Z2) - (Y - Z2);
1963 X2 = X8 / (Y - Half);
1964 X2 = X2 - Half * X2 * X2;
1965 SqEr = (Y2 + Half) + (Half - X2);
1966 if (SqEr < MinSqEr) MinSqEr = SqEr;
1968 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1978 if (Z == Zero && Q <= Zero)
1979 printf("WARNING: computing\n");
1980 else BadCond(Defect, "computing\n");
1981 printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
1982 printf("\tyielded %.17e;\n", Y);
1983 printf("\twhich compared unequal to correct %.17e ;\n",
1985 printf("\t\tthey differ by %.17e .\n", Y - X);
1987 N = N + 1; /* ... count discrepancies. */
2004 /* PrintIfNPositive */
2008 if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
2017 printf("Since comparison denies Z = 0, evaluating ");
2018 printf("(Z + Z) / Z should be safe.\n");
2020 if (setjmp(ovfl_buf)) goto very_serious;
2022 printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
2024 if (FABS(Q9 - Two) < Radix * U2) {
2025 printf("This is O.K., provided Over/Underflow");
2026 printf(" has NOT just been signaled.\n");
2029 if ((Q9 < One) || (Q9 > Two)) {
2032 ErrCnt [Serious] = ErrCnt [Serious] + 1;
2033 printf("This is a VERY SERIOUS DEFECT!\n");
2037 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2038 printf("This is a DEFECT!\n");
2047 if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
2052 BadCond(Defect, "What prints as Z = ");
2053 printf("%.17e\n\tcompares different from ", Z);
2054 if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
2055 if (! ((Z == Random2)
2056 || (Random2 == Random1)))
2057 printf("1 * Z == %g\n", Random2);
2058 if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
2059 if (Random2 != Random1) {
2060 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2061 BadCond(Defect, "Multiplication does not commute!\n");
2062 printf("\tComparison alleges that 1 * Z = %.17e\n",
2064 printf("\tdiffers from Z * 1 = %.17e\n", Random1);
2074 printf("%s test appears to be inconsistent...\n", s);
2075 printf(" PLEASE NOTIFY KARPINKSI!\n");
2084 { while(*s) printf("%s\n", *s++); }
2088 static char *instr[] = {
2089 "Lest this program stop prematurely, i.e. before displaying\n",
2090 " `END OF TEST',\n",
2091 "try to persuade the computer NOT to terminate execution when an",
2092 "error like Over/Underflow or Division by Zero occurs, but rather",
2093 "to persevere with a surrogate value after, perhaps, displaying some",
2094 "warning. If persuasion avails naught, don't despair but run this",
2095 "program anyway to see how many milestones it passes, and then",
2096 "amend it to make further progress.\n",
2097 "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
2107 static char *head[] = {
2108 "Users are invited to help debug and augment this program so it will",
2109 "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
2110 "Please send suggestions and interesting results to",
2111 "\tRichard Karpinski",
2112 "\tComputer Center U-76",
2113 "\tUniversity of California",
2114 "\tSan Francisco, CA 94143-0704, USA\n",
2115 "In doing so, please include the following information:",
2117 "\tPrecision:\tsingle;",
2119 "\tPrecision:\tdouble;",
2121 "\tVersion:\t10 February 1989;",
2124 "\tOptimization level:\n",
2125 "\tOther relevant compiler options:",
2131 /* Characteristics */
2135 static char *chars[] = {
2136 "Running this program should reveal these characteristics:",
2137 " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
2138 " Precision = number of significant digits carried.",
2139 " U2 = Radix/Radix^Precision = One Ulp",
2140 "\t(OneUlpnit in the Last Place) of 1.000xxx .",
2141 " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
2142 " Adequacy of guard digits for Mult., Div. and Subt.",
2143 " Whether arithmetic is chopped, correctly rounded, or something else",
2144 "\tfor Mult., Div., Add/Subt. and Sqrt.",
2145 " Whether a Sticky Bit used correctly for rounding.",
2146 " UnderflowThreshold = an underflow threshold.",
2147 " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
2148 " V = an overflow threshold, roughly.",
2149 " V0 tells, roughly, whether Infinity is represented.",
2150 " Comparisions are checked for consistency with subtraction",
2151 "\tand for contamination with pseudo-zeros.",
2152 " Sqrt is tested. Y^X is not tested.",
2153 " Extra-precise subexpressions are revealed but NOT YET tested.",
2154 " Decimal-Binary conversion is NOT YET tested for accuracy.",
2163 /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
2164 with further massaging by David M. Gay. */
2166 static char *hist[] = {
2167 "The program attempts to discriminate among",
2168 " FLAWs, like lack of a sticky bit,",
2169 " Serious DEFECTs, like lack of a guard digit, and",
2170 " FAILUREs, like 2+2 == 5 .",
2171 "Failures may confound subsequent diagnoses.\n",
2172 "The diagnostic capabilities of this program go beyond an earlier",
2173 "program called `MACHAR', which can be found at the end of the",
2174 "book `Software Manual for the Elementary Functions' (1980) by",
2175 "W. J. Cody and W. Waite. Although both programs try to discover",
2176 "the Radix, Precision and range (over/underflow thresholds)",
2177 "of the arithmetic, this program tries to cope with a wider variety",
2178 "of pathologies, and to say how well the arithmetic is implemented.",
2179 "\nThe program is based upon a conventional radix representation for",
2180 "floating-point numbers, but also allows logarithmic encoding",
2181 "as used by certain early WANG machines.\n",
2182 "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
2183 "see source comments for more history.",
2190 pow(x, y) /* return x ^ y (exponentiation) */
2193 extern double exp(), frexp(), ldexp(), log(), modf();
2196 int ex, ey = 0, flip = 0;
2200 if ((y < -1100. || y > 1100.) && x != -1.) return exp(y * log(x));
2202 if (y < 0.) { y = -y; flip = 1; }
2204 if (y) xy = exp(y * log(x));
2206 /* next several lines assume >= 32 bit integers */
2208 if (i = ye) for(;;) {
2209 if (i & 1) { xy *= x; ey += ex; }
2210 if (!(i >>= 1)) break;
2213 if (x < .5) { x *= 2.; ex -= 1; }
2215 if (flip) { xy = 1. / xy; ey = -ey; }
2216 return ldexp(xy, ey);
2219 #endif /* NO_FLOATS */