2 !!DESCRIPTION!! a wellknown floating point test
3 !!ORIGIN!! LCC 4.1 Testsuite
4 !!LICENCE!! own, freely distributeable for non-profit. read CPYRIGHT.LCC
11 printf("NO_FLOATS\n\r");
18 /* A C version of Kahan's Floating Point Test "Paranoia"
20 Thos Sumner, UCSF, Feb. 1985
21 David Gay, BTL, Jan. 1986
23 This is a rewrite from the Pascal version by
25 B. A. Wichmann, 18 Jan. 1985
27 (and does NOT exhibit good C programming style).
29 (C) Apr 19 1983 in BASIC version by:
30 Professor W. M. Kahan,
32 Electrical Engineering & Computer Science Dept.
33 University of California
34 Berkeley, California 94720
37 converted to Pascal by:
39 National Physical Laboratory
46 David M. Gay and Thos Sumner
47 AT&T Bell Labs Computer Center, Rm. U-76
48 600 Mountain Avenue University of California
49 Murray Hill, NJ 07974 San Francisco, CA 94143
52 with simultaneous corrections to the Pascal source (reflected
53 in the Pascal source available over netlib).
54 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
56 Reports of results on various systems from all the versions
57 of Paranoia are being collected by Richard Karpinski at the
58 same address as Thos Sumner. This includes sample outputs,
59 bug reports, and criticisms.
61 You may copy this program freely if you acknowledge its source.
62 Comments on the Pascal version to NPL, please.
64 The C version catches signals from floating-point exceptions.
65 If signal(SIGFPE,...) is unavailable in your environment, you may
66 #define NOSIGNAL to comment out the invocations of signal.
68 This source file is too big for some C compilers, but may be split
69 into pieces. Comments containing "SPLIT" suggest convenient places
70 for this splitting. At the end of these comments is an "ed script"
71 (for the UNIX(tm) editor ed) that will do this splitting.
73 By #defining Single when you compile this source, you may obtain
74 a single-precision C version of Paranoia.
76 The following is from the introductory commentary from Wichmann's work:
78 The BASIC program of Kahan is written in Microsoft BASIC using many
79 facilities which have no exact analogy in Pascal. The Pascal
80 version below cannot therefore be exactly the same. Rather than be
81 a minimal transcription of the BASIC program, the Pascal coding
82 follows the conventional style of block-structured languages. Hence
83 the Pascal version could be useful in producing versions in other
86 Rather than use identifiers of minimal length (which therefore have
87 little mnemonic significance), the Pascal version uses meaningful
88 identifiers as follows [Note: A few changes have been made for C]:
90 BASIC C BASIC C BASIC C
93 A1 AInverse J0 NoErrors T
94 B Radix [Failure] T0 Underflow
95 B1 BInverse J1 NoErrors T2 ThirtyTwo
96 B2 RadixD2 [SeriousDefect] T5 OneAndHalf
97 B9 BMinusU2 J2 NoErrors T7 TwentySeven
98 C [Defect] T8 TwoForty
99 C1 CInverse J3 NoErrors U OneUlp
100 D [Flaw] U0 UnderflowThreshold
111 F1 MinusOne O5 Five X9 Random1
117 G2 GDiv Q9 Z0 PseudoZero
120 H1 HInverse R2 RDiv Z9
127 All the variables in BASIC are true variables and in consequence,
128 the program is more difficult to follow since the "constants" must
129 be determined (the glossary is very helpful). The Pascal version
130 uses Real constants, but checks are added to ensure that the values
131 are correctly converted by the compiler.
133 The major textual change to the Pascal version apart from the
134 identifiersis that named procedures are used, inserting parameters
135 wherehelpful. New procedures are also introduced. The
136 correspondence is as follows:
142 170- 250 Instructions
144 480- 670 Characteristics
148 4040-4080 DoesYequalX
149 4090-4110 PrintIfNPositive
150 4640-4850 TestPartialUnderflow
152 =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
154 Below is an "ed script" that splits para.c into 10 files
155 of the form part[1-8].c, subs.c, and msgs.c, plus a header
156 file, paranoia.h, that these files require.
213 /Computed constants/,$d
215 1,$s/^FLOAT/extern &/
218 /^Guard/,/^Round/s/^/extern /
219 /^jmp_buf/s/^/extern /
220 /^Sig_type/s/^/extern /
222 extern void sigfpe();/
234 extern double fabs(), floor(), log(), pow(), sqrt();
238 #define FABS(x) (float)fabs((double)(x))
239 #define FLOOR(x) (float)floor((double)(x))
240 #define LOG(x) (float)log((double)(x))
241 #define POW(x,y) (float)pow((double)(x),(double)(y))
242 #define SQRT(x) (float)sqrt((double)(x))
245 #define FABS(x) fabs(x)
246 #define FLOOR(x) floor(x)
247 #define LOG(x) log(x)
248 #define POW(x,y) pow(x,y)
249 #define SQRT(x) sqrt(x)
253 typedef void (*Sig_type)();
258 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
259 FLOAT Sign(), Random();
261 /*Small floating point constants.*/
271 FLOAT TwentySeven = 27.0;
272 FLOAT ThirtyTwo = 32.0;
273 FLOAT TwoForty = 240.0;
274 FLOAT MinusOne = -1.0;
275 FLOAT OneAndHalf = 1.5;
276 /*Integer constants*/
277 int NoTrials = 20; /*Number of tests for commutativity. */
281 /* Definitions for declared types
283 Rounding == (Chopped, Rounded, Other);
284 Message == packed array [1..40] of char;
285 Class == (Flaw, Defect, Serious, Failure);
296 typedef int Guard, Rounding, Class;
297 typedef char Message;
299 /* Declarations of Variables */
305 FLOAT E0, E1, Exp2, E3, MinSqEr;
306 FLOAT SqEr, MaxSqEr, E9;
316 FLOAT T, Underflow, S;
317 FLOAT OneUlp, UfThold, U1, U2;
320 FLOAT X, X1, X2, X8, Random1;
321 FLOAT Y, Y1, Y2, Random2;
322 FLOAT Z, PseudoZero, Z1, Z2, Z9;
328 Guard GMult, GDiv, GAddSub;
329 Rounding RMult, RDiv, RAddSub, RSqrt;
330 int Break, Done, NotMonot, Monot, Anomaly, IEEE,
332 /* Computed constants. */
333 /*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
334 /*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
336 /* floating point exception receiver */
341 printf("\n* * * FLOATING-POINT ERROR * * *\n");
345 signal(SIGFPE, sigsave);
348 longjmp(ovfl_buf, 1);
357 ieee_flags("set", "precision", "double", &out);
359 /* First two assignments use integer right-hand sides. */
367 Nine = Three * Three;
368 TwentySeven = Nine * Three;
369 ThirtyTwo = Four * Eight;
370 TwoForty = Four * Five * Three * Four;
373 OneAndHalf = One + Half;
379 /*=============================================*/
381 /*=============================================*/
383 signal(SIGFPE, sigfpe);
393 /*=============================================*/
395 /*=============================================*/
396 printf("Program is now RUNNING tests on small integers:\n");
398 TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
399 && (One > Zero) && (One + One == Two),
400 "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
403 ErrCnt[Failure] = ErrCnt[Failure] + 1;
404 printf("Comparison alleges that -0.0 is Non-zero!\n");
409 TstCond (Failure, (Three == Two + One) && (Four == Three + One)
410 && (Four + Two * (- Two) == Zero)
411 && (Four - Three - One == Zero),
412 "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
413 TstCond (Failure, (MinusOne == (0 - One))
414 && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
415 && (MinusOne + FABS(One) == Zero)
416 && (MinusOne + MinusOne * MinusOne == Zero),
417 "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
418 TstCond (Failure, Half + MinusOne + Half == Zero,
419 "1/2 + (-1) + 1/2 != 0");
420 /*=============================================*/
430 #include "paranoia.h"
434 /*=============================================*/
435 TstCond (Failure, (Nine == Three * Three)
436 && (TwentySeven == Nine * Three) && (Eight == Four + Four)
437 && (ThirtyTwo == Eight * Four)
438 && (ThirtyTwo - TwentySeven - Four - One == Zero),
439 "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
440 TstCond (Failure, (Five == Four + One) &&
441 (TwoForty == Four * Five * Three * Four)
442 && (TwoForty / Three - Four * Four * Five == Zero)
443 && ( TwoForty / Four - Five * Three * Four == Zero)
444 && ( TwoForty / Five - Four * Three * Four == Zero),
445 "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
446 if (ErrCnt[Failure] == 0) {
447 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
450 printf("Searching for Radix and Precision.\n");
457 } while (MinusOne + FABS(Y) < Zero);
458 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
465 } while ( Radix == Zero);
466 if (Radix < Two) Radix = One;
467 printf("Radix = %f .\n", Radix);
471 Precision = Precision + One;
474 } while ((Y - W) == One);
476 /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
480 printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
481 printf("Recalculating radix and precision\n ");
496 /*... now X = (unknown no.) ulps of 1+...*/
499 Y = Half * U2 + ThirtyTwo * U2 * U2;
502 } while ( ! ((U2 <= X) || (X <= Zero)));
504 /*... now U2 == 1 ulp of 1 + ... */
512 /*... now X == (unknown no.) ulps of 1 -... */
515 Y = Half * U1 + ThirtyTwo * U1 * U1;
520 } while ( ! ((U1 <= X) || (X <= Zero)));
521 /*... now U1 == 1 ulp of 1 - ... */
522 if (U1 == E1) printf("confirms closest relative separation U1 .\n");
523 else printf("gets better closest relative separation U1 = %.7e .\n", U1);
525 F9 = (Half - U1) + Half;
526 Radix = FLOOR(0.01 + U2 / U1);
527 if (Radix == E0) printf("Radix confirmed.\n");
528 else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
529 TstCond (Defect, Radix <= Eight + Eight,
530 "Radix is too big: roundoff problems");
531 TstCond (Flaw, (Radix == Two) || (Radix == 10)
532 || (Radix == One), "Radix is not as good as 2 or 10");
533 /*=============================================*/
535 /*=============================================*/
536 TstCond (Failure, F9 - Half < Half,
537 "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
542 TstCond (Failure, (X != One)
543 || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
546 /*=============================================*/
548 /*=============================================*/
549 /*... BMinusU2 = nextafter(Radix, 0) */
550 BMinusU2 = Radix - One;
551 BMinusU2 = (BMinusU2 - U2) + One;
552 /* Purify Integers */
554 X = - TwoForty * LOG(U1) / LOG(Radix);
556 if (FABS(X - Y) * Four < One) X = Y;
557 Precision = X / TwoForty;
558 Y = FLOOR(Half + Precision);
559 if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
561 if ((Precision != FLOOR(Precision)) || (Radix == One)) {
562 printf("Precision cannot be characterized by an Integer number\n");
563 printf("of significant digits but, by itself, this is a minor flaw.\n");
566 printf("logarithmic encoding has precision characterized solely by U1.\n");
567 else printf("The number of significant digits of the Radix is %f .\n",
569 TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
570 "Precision worse than 5 decimal figures ");
571 /*=============================================*/
573 /*=============================================*/
574 /* Test for extra-precise subepressions */
575 X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
578 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
579 } while ( ! ((Z2 <= X) || (X <= Zero)));
580 X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
583 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
584 + One / Two)) + One / Two;
585 } while ( ! ((Z1 <= Z) || (Z <= Zero)));
589 Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
591 } while ( ! ((Y1 <= Y) || (Y <= Zero)));
593 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
594 } while ( ! ((X1 <= X) || (X <= Zero)));
595 if ((X1 != Y1) || (X1 != Z1)) {
596 BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
597 printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
598 printf("are symptoms of inconsistencies introduced\n");
599 printf("by extra-precise evaluation of arithmetic subexpressions.\n");
600 notify("Possibly some part of this");
601 if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
602 "That feature is not tested further by this program.\n") ;
605 if ((Z1 != U1) || (Z2 != U2)) {
606 if ((Z1 >= U1) || (Z2 >= U2)) {
607 BadCond(Failure, "");
609 printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
610 printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
613 if ((Z1 <= Zero) || (Z2 <= Zero)) {
614 printf("Because of unusual Radix = %f", Radix);
615 printf(", or exact rational arithmetic a result\n");
616 printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
617 notify("of an\nextra-precision");
619 if (Z1 != Z2 || Z1 > Zero) {
624 printf("Some subexpressions appear to be calculated extra\n");
625 printf("precisely with about %g extra B-digits, i.e.\n",
627 printf("roughly %g extra significant decimals.\n",
630 printf("That feature is not tested further by this program.\n");
635 /*=============================================*/
638 #include "paranoia.h"
642 /*=============================================*/
644 X = W / (Radix * Radix);
649 TstCond (Failure, X == U2,
650 "Subtraction is not normalized X=Y,X+Z != Y+Z!");
652 "Subtraction appears to be normalized, as it should be.");
654 printf("\nChecking for guard digit in *, /, and -.\n");
667 X = X * (Radix - One);
668 T = T * (Radix - One);
669 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
672 TstCond (Serious, False,
673 "* lacks a Guard Digit, so 1*X != X");
677 Y = FABS((X + Z) - X * X) - U2;
679 Z = FABS((X - U2) - X * X) - U1;
680 TstCond (Failure, (Y <= Zero)
681 && (Z <= Zero), "* gets too many final digits wrong.\n");
689 T = Nine / TwentySeven;
691 TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
692 "Division lacks a Guard Digit, so error can exceed 1 ulp\nor 1/3 and 3/9 and 9/27 may disagree");
699 if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
702 TstCond (Serious, False,
703 "Division lacks a Guard Digit, so X/1 != X");
705 X = One / (One + U2);
707 TstCond (Serious, Y < Zero,
708 "Computed value of 1/1.000..1 >= 1");
710 Y = One + Radix * U2;
714 StickyBit = T / Radix;
717 TstCond (Failure, X == Zero && Y == Zero,
718 "* and/or / gets too many last digits wrong");
723 Z = Radix - BMinusU2;
725 if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
728 TstCond (Serious, False,
729 "- lacks Guard Digit, so cancellation is obscured");
731 if (F9 != One && F9 - One >= Zero) {
732 BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
733 printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
734 printf(" such precautions against division by zero as\n");
735 printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
737 if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
738 " *, /, and - appear to have guard digits, as they should.\n");
739 /*=============================================*/
741 /*=============================================*/
743 printf("Checking rounding on multiply, divide and add/subtract.\n");
747 RadixD2 = Radix / Two;
754 AInvrse = AInvrse / A1;
755 } while ( ! (FLOOR(AInvrse) != AInvrse));
756 Done = (X == One) || (A1 > Three);
757 if (! Done) A1 = Nine + One;
759 if (X == One) A1 = Radix;
766 TstCond (Failure, Z == Half,
767 "X * (1/X) differs from 1");
784 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
785 X = (OneAndHalf + U2) * Y2;
786 Y = OneAndHalf - U2 - U2;
787 Z = OneAndHalf + U2 + U2;
788 T = (OneAndHalf - U2) * Y1;
793 Y = (U2 - Y) + StickyBit;
794 Z = S - (Z + U2 + U2);
795 StickyBit = (Y2 + U2) * Y1;
797 StickyBit = StickyBit - Y2;
799 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
800 && ( StickyBit == Zero) && (Y1 == Half)) {
802 printf("Multiplication appears to round correctly.\n");
804 else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
805 && (T < Zero) && (StickyBit + U2 == Zero)
808 printf("Multiplication appears to chop.\n");
810 else printf("* is neither chopped nor correctly rounded.\n");
811 if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
813 else printf("* is neither chopped nor correctly rounded.\n");
814 /*=============================================*/
816 /*=============================================*/
819 Z = OneAndHalf + U2 + U2;
821 T = OneAndHalf - U2 - U2;
827 Z = Z - (OneAndHalf + U2);
828 T = (U2 - OneAndHalf) + T;
829 if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
843 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
844 && (Y2 == Zero) && (Y2 == Zero)
845 && (Y1 - Half == F9 - Half )) {
847 printf("Division appears to round correctly.\n");
848 if (GDiv == No) notify("Division");
850 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
851 && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
853 printf("Division appears to chop.\n");
856 if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
857 BInvrse = One / Radix;
858 TstCond (Failure, (BInvrse * Radix - Half == Half),
859 "Radix * ( 1 / Radix ) differs from 1");
860 /*=============================================*/
863 #include "paranoia.h"
867 /*=============================================*/
868 TstCond (Failure, ((F9 + U1) - Half == Half)
869 && ((BMinusU2 + U2 ) - One == Radix - One),
870 "Incomplete carry-propagation in Addition");
872 Y = One + U2 * (One - U2);
876 if ((X == Zero) && (Y == Zero)) {
878 printf("Add/Subtract appears to be chopped.\n");
880 if (GAddSub == Yes) {
881 X = (Half + U2) * U2;
882 Y = (Half - U2) * U2;
887 if ((X == Zero) && (Y == Zero)) {
888 X = (Half + U2) * U1;
889 Y = (Half - U2) * U1;
894 if ((X == Zero) && (Y == Zero)) {
896 printf("Addition/Subtraction appears to round correctly.\n");
897 if (GAddSub == No) notify("Add/Subtract");
899 else printf("Addition/Subtraction neither rounds nor chops.\n");
901 else printf("Addition/Subtraction neither rounds nor chops.\n");
903 else printf("Addition/Subtraction neither rounds nor chops.\n");
905 X = One + Half * (One + Half);
906 Y = (One + U2) * Half;
910 if (StickyBit != Zero) {
912 BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
915 if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
916 && (RMult == Rounded) && (RDiv == Rounded)
917 && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
918 printf("Checking for sticky bit.\n");
919 X = (Half + U1) * U2;
923 if ((Z - One <= Zero) && (T - One >= U2)) {
926 if ((Z - T >= U2) && (Y - T == Zero)) {
927 X = (Half + U1) * U1;
931 if ((Z - One == Zero) && (T - F9 == Zero)) {
932 Z = (Half - U1) * U1;
935 if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
936 Z = (One + U2) * OneAndHalf;
937 T = (OneAndHalf + U2) - Z + U2;
938 X = One + Half / Radix;
939 Y = One + Radix * U2;
941 if (T == Zero && X + Radix * U2 - Z == Zero) {
945 if ((Y - One == Zero)) StickyBit = S;
954 if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
955 else printf("Sticky bit used incorrectly or not at all.\n");
956 TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
957 RMult == Other || RDiv == Other || RAddSub == Other),
958 "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");
959 /*=============================================*/
961 /*=============================================*/
963 printf("Does Multiplication commute? ");
964 printf("Testing on %d random pairs.\n", NoTrials);
975 } while ( ! ((I > NoTrials) || (Z9 != Zero)));
977 Random1 = One + Half / Three;
978 Random2 = (U2 + U1) + One;
979 Z = Random1 * Random2;
980 Y = Random2 * Random1;
981 Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
982 Three) * ((U2 + U1) + One);
984 if (! ((I == NoTrials) || (Z9 == Zero)))
985 BadCond(Defect, "X * Y == Y * X trial fails.\n");
986 else printf(" No failures found in %d integer pairs.\n", NoTrials);
987 /*=============================================*/
989 /*=============================================*/
990 printf("\nRunning test of square root(x).\n");
991 TstCond (Failure, (Zero == SQRT(Zero))
992 && (- Zero == SQRT(- Zero))
993 && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
1001 OneUlp = BInvrse * U1;
1006 if (J != Zero) Pause();
1007 printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
1011 if ((Radix != One)) do {
1014 } while ( ! ((Y - X >= NoTrials)));
1017 while (I <= NoTrials) {
1020 if (J > Zero) break;
1023 printf("Test for sqrt monotonicity.\n");
1027 Z = Radix + Radix * U2;
1030 while ( ! (NotMonot || Monot)) {
1035 if ((X > Q) || (Q > Z)) NotMonot = True;
1037 Q = FLOOR(Q + Half);
1038 if ((I > 0) || (Radix == Q * Q)) Monot = True;
1040 if (I > 1) Monot = True;
1054 if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
1056 BadCond(Defect, "");
1057 printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
1059 /*=============================================*/
1062 #include "paranoia.h"
1066 /*=============================================*/
1067 MinSqEr = MinSqEr + Half;
1068 MaxSqEr = MaxSqEr - Half;
1069 Y = (SQRT(One + U2) - One) / U2;
1070 SqEr = (Y - One) + U2 / Eight;
1071 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1072 SqEr = Y + U2 / Eight;
1073 if (SqEr < MinSqEr) MinSqEr = SqEr;
1074 Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
1075 SqEr = Y + U1 / Eight;
1076 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1077 SqEr = (Y + One) + U1 / Eight;
1078 if (SqEr < MinSqEr) MinSqEr = SqEr;
1081 for( Indx = 1; Indx <= 3; ++Indx) {
1082 Y = SQRT((X + U1 + X) + F9);
1083 Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
1084 Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
1085 SqEr = (Y + Half) + Z;
1086 if (SqEr < MinSqEr) MinSqEr = SqEr;
1087 SqEr = (Y - Half) + Z;
1088 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1089 if (((Indx == 1) || (Indx == 3)))
1090 X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
1096 /*=============================================*/
1098 /*=============================================*/
1101 RSqrt = Other; /* ~dgh */
1103 printf("Testing whether sqrt is rounded or chopped.\n");
1104 D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
1105 /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
1108 if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
1123 X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
1126 } while ( ! (X1 <= Zero));
1135 if (Y2 >= FourD) Y2 = Y2 - FourD;
1136 } while ( ! (Y >= D));
1138 Q = (X8 + Z * Z) / FourD;
1140 if (Q != FLOOR(Q)) Anomaly = True;
1145 X = X - FLOOR(X / Radix) * Radix;
1150 } while ( ! (Break || (Z1 <= Zero)));
1151 if ((Z1 <= Zero) && (! Break)) Anomaly = True;
1153 if (Z1 > RadixD2) Z1 = Z1 - Radix;
1156 } while ( ! (U2 * D >= F9));
1157 if (D * Radix - D != W - D) Anomaly = True;
1161 Y = D + (One + Z) * Half;
1164 Y = D + (One - Z) * Half + D;
1169 if (D - Z2 != W - Z2) Anomaly = True;
1171 Y = (D - Z2) + (Z2 + (One - Z) * Half);
1172 X = (D - Z2) + (Z2 - Z + Q);
1174 Y = (One + Z) * Half;
1177 if (I == 0) Anomaly = True;
1183 if ((I == 0) || Anomaly) {
1184 BadCond(Failure, "Anomalous arithmetic with Integer < ");
1185 printf("Radix^Precision = %.7e\n", W);
1186 printf(" fails test whether sqrt rounds or chops.\n");
1191 if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
1193 printf("Square root appears to be correctly rounded.\n");
1196 if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
1197 || (MinSqEr + Radix < Half)) SqRWrng = True;
1200 printf("Square root appears to be chopped.\n");
1205 printf("Square root is neither chopped nor correctly rounded.\n");
1206 printf("Observed errors run from %.7e ", MinSqEr - Half);
1207 printf("to %.7e ulps.\n", Half + MaxSqEr);
1208 TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
1209 "sqrt gets too many last digits wrong");
1211 /*=============================================*/
1213 /*=============================================*/
1215 printf("Testing powers Z^i for small Integers Z and i.\n");
1217 /* ... test powers of zero. */
1229 if (Z == MinusOne) Break = True;
1234 /* .. if(-1)^N is invalid, replace MinusOne by One. */
1242 M = FLOOR(Two * LOG(W) / LOG(A1));
1248 if (Z == AInvrse) Break = True;
1250 } while ( ! (Break));
1251 /*=============================================*/
1253 /*=============================================*/
1254 /* Powers of Radix have been tested, */
1255 /* next try a few primes */
1264 } while ( Three * FLOOR(Z / Three) == Z );
1265 } while ( Z < Eight * Three );
1267 printf("Errors like this may invalidate financial calculations\n");
1268 printf("\tinvolving interest rates.\n");
1272 if (N == 0) printf("... no discrepancis found.\n");
1275 /*=============================================*/
1278 #include "paranoia.h"
1282 /*=============================================*/
1283 printf("Seeking Underflow thresholds UfThold and E0.\n");
1285 if (Precision != FLOOR(Precision)) {
1291 } while ( X > Zero);
1295 /* ... D is power of 1/Radix < 1. */
1300 } while ((Y > Z) && (Z + Z > Z));
1307 } while ((Y > Z) && (Z + Z > Z));
1308 if (Radix < Two) HInvrse = Two;
1309 else HInvrse = Radix;
1311 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1315 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1320 } while ((E0 > Z) && (Z + Z > Z));
1332 BadCond(Failure, "multiplication gets too many last digits wrong.\n");
1341 PseudoZero = Underflow * H;
1345 Underflow = PseudoZero;
1346 if (E1 + E1 <= E1) {
1347 Y2 = Underflow * HInvrse;
1350 if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
1352 PseudoZero = PseudoZero * H;
1353 } while ((Underflow > PseudoZero)
1354 && (PseudoZero + PseudoZero > PseudoZero));
1356 /* Comment line 4530 .. 4560 */
1357 if (PseudoZero != Zero) {
1360 /* ... Test PseudoZero for "phoney- zero" violates */
1361 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1363 if (PseudoZero <= Zero) {
1364 BadCond(Failure, "Positive expressions can underflow to an\n");
1365 printf("allegedly negative value\n");
1366 printf("PseudoZero that prints out as: %g .\n", PseudoZero);
1369 printf("But -PseudoZero, which should be\n");
1370 printf("positive, isn't; it prints out as %g .\n", X);
1374 BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
1375 printf("value PseudoZero that prints out as %g .\n", PseudoZero);
1379 /*=============================================*/
1381 /*=============================================*/
1382 if (CInvrse * Y > CInvrse * Y1) {
1386 if (! ((E1 == Zero) || (E1 == E0))) {
1387 BadCond(Defect, "");
1389 printf("Products underflow at a higher");
1390 printf(" threshold than differences.\n");
1391 if (PseudoZero == Zero)
1395 printf("Difference underflows at a higher");
1396 printf(" threshold than products.\n");
1399 printf("Smallest strictly positive number found is E0 = %g .\n", E0);
1403 if (N == 1) Underflow = Y;
1405 if (E1 == Zero) I = 3;
1406 if (UfThold == Zero) I = I - 2;
1410 UfThold = Underflow;
1411 if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
1413 BadCond(Failure, "Either accuracy deteriorates as numbers\n");
1414 printf("approach a threshold = %.17e\n", UfThold);;
1415 printf(" coming down from %.17e\n", C);
1416 printf(" or else multiplication gets too many last digits wrong.\n");
1422 BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
1423 printf("Q == Y while denying that |Q - Y| == 0; these values\n");
1424 printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
1425 printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
1434 if ((Q == UfThold) && (E1 == E0)
1435 && (FABS( UfThold - E1 / E9) <= E1)) {
1437 printf("Underflow is gradual; it incurs Absolute Error =\n");
1438 printf("(roundoff in UfThold) < E0.\n");
1440 Y = Y * (OneAndHalf + U2);
1441 X = CInvrse * (One + U2);
1449 if (setjmp(ovfl_buf)) {
1450 printf("Underflow / UfThold failed!\n");
1453 else R = SQRT(Underflow / UfThold);
1457 X = Z * (One + R * H * (One + H));
1461 X = Z * (One + H * H * (One + H));
1463 if (! ((X == Z) || (X - Z != Zero))) {
1465 printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
1467 printf("yet X - Z yields %.17e .\n", Z9);
1468 printf(" Should this NOT signal Underflow, ");
1469 printf("this is a SERIOUS DEFECT\nthat causes ");
1470 printf("confusion when innocent statements like\n");;
1471 printf(" if (X == Z) ... else");
1472 printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
1473 printf("encounter Division by Zero although actually\n");
1475 if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
1476 else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
1480 printf("The Underflow threshold is %.17e, %s\n", UfThold,
1482 printf("calculation may suffer larger Relative error than ");
1483 printf("merely roundoff.\n");
1487 if (Y2 <= UfThold) {
1489 BadCond(Defect, "");
1493 BadCond(Serious, "");
1496 printf("Range is too narrow; U1^%d Underflows.\n", I);
1498 /*=============================================*/
1501 #include "paranoia.h"
1505 /*=============================================*/
1506 Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
1508 printf("Since underflow occurs below the threshold\n");
1509 printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
1510 printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
1511 V9 = POW(HInvrse, Y2);
1512 printf("actually calculating yields: %.17e .\n", V9);
1513 if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
1514 BadCond(Serious, "this is not between 0 and underflow\n");
1515 printf(" threshold = %.17e .\n", UfThold);
1517 else if (! (V9 > UfThold * (One + E9)))
1518 printf("This computed value is O.K.\n");
1520 BadCond(Defect, "this is not between 0 and underflow\n");
1521 printf(" threshold = %.17e .\n", UfThold);
1523 /*=============================================*/
1525 /*=============================================*/
1527 /* ...calculate Exp2 == exp(2) == 7.389056099... */
1541 Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
1546 printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
1550 Z = (X + One) / (Z - (One - BInvrse));
1551 Q = POW(X, Z) - Exp2;
1552 if (FABS(Q) > TwoForty * U2) {
1554 V9 = (X - BInvrse) - (One - BInvrse);
1555 BadCond(Defect, "Calculated");
1556 printf(" %.17e for\n", POW(X,Z));
1557 printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
1558 printf("\tdiffers from correct value by %.17e .\n", Q);
1559 printf("\tThis much error may spoil financial\n");
1560 printf("\tcalculations involving tiny interest rates.\n");
1564 Z = (Y - X) * Two + Y;
1567 Z = One + (X - F9)*(X - F9);
1568 if (Z > One && I < NoTrials) I++;
1572 printf("Accuracy seems adequate.\n");
1584 /*=============================================*/
1586 /*=============================================*/
1587 printf("Testing powers Z^Q at four nearly extreme values.\n");
1590 Q = FLOOR(Half - LOG(C) / LOG(A1));
1600 if (Z < One) Break = True;
1602 } while ( ! (Break));
1604 if (N == 0) printf(" ... no discrepancies found.\n");
1607 /*=============================================*/
1609 /*=============================================*/
1611 printf("Searching for Overflow threshold:\n");
1612 printf("This may generate an error.\n");
1616 if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
1626 printf("Can `Z = -Y' overflow?\n");
1627 printf("Trying it on Y = %.17e .\n", Y);
1630 if (V - Y == V + V0) printf("Seems O.K.\n");
1633 BadCond(Flaw, "-(-Y) differs from Y.\n");
1636 BadCond(Serious, "");
1637 printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
1640 Y = V * (HInvrse * U2 - HInvrse);
1641 Z = Y + ((One - HInvrse) * U2) * V;
1644 if (V0 - V < V0) V = V0;
1647 V = Y * (HInvrse * U2 - HInvrse);
1648 V = V + ((One - HInvrse) * U2) * Y;
1650 printf("Overflow threshold is V = %.17e .\n", V);
1651 if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
1652 else printf("There is no saturation value because the system traps on overflow.\n");
1654 printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
1656 printf(" nor for V / 1 = %.17e .\n", V9);
1657 printf("Any overflow signal separating this * from the one\n");
1658 printf("above is a DEFECT.\n");
1659 /*=============================================*/
1661 /*=============================================*/
1662 if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
1663 BadCond(Failure, "Comparisons involving ");
1664 printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
1667 /*=============================================*/
1669 /*=============================================*/
1671 for(Indx = 1; Indx <= 3; ++Indx) {
1673 case 1: Z = UfThold; break;
1674 case 2: Z = E0; break;
1675 case 3: Z = PseudoZero; break;
1680 if (Y / (One - Radix * E9) < Z
1681 || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
1682 if (V9 > U1) BadCond(Serious, "");
1683 else BadCond(Defect, "");
1684 printf("Comparison alleges that what prints as Z = %.17e\n", Z);
1685 printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
1689 /*=============================================*/
1691 /*=============================================*/
1692 for(Indx = 1; Indx <= 2; ++Indx) {
1693 if (Indx == 1) Z = V;
1696 X = (One - Radix * E9) * V9;
1698 if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
1700 if (X < W) BadCond(Serious, "");
1701 else BadCond(Defect, "");
1702 printf("Comparison alleges that Z = %17e\n", Z);
1703 printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
1706 /*=============================================*/
1709 #include "paranoia.h"
1713 /*=============================================*/
1717 if (X*Y < One || X > Y) {
1718 if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
1719 else BadCond(Flaw, "");
1721 printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
1722 X, "is too far from 1.\n");
1724 /*=============================================*/
1726 /*=============================================*/
1727 for (Indx = 1; Indx <= 5; ++Indx) {
1730 case 2: X = One + U2; break;
1731 case 3: X = V; break;
1732 case 4: X = UfThold; break;
1737 if (setjmp(ovfl_buf))
1738 printf(" X / X traps when X = %g\n", X);
1740 V9 = (Y / X - Half) - Half;
1741 if (V9 == Zero) continue;
1742 if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
1743 else BadCond(Serious, "");
1744 printf(" X / X differs from 1 when X = %.17e\n", X);
1745 printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
1749 /*=============================================*/
1751 /*=============================================*/
1754 printf("What message and/or values does Division by Zero produce?\n") ;
1756 printf("This can interupt your program. You can ");
1757 printf("skip this part if you wish.\n");
1758 printf("Do you wish to compute 1 / 0? ");
1760 read (KEYBOARD, ch, 8);
1761 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1764 printf(" Trying to compute 1 / 0 produces ...");
1765 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
1769 else printf("O.K.\n");
1770 printf("\nDo you wish to compute 0 / 0? ");
1772 read (KEYBOARD, ch, 80);
1773 if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1776 printf("\n Trying to compute 0 / 0 produces ...");
1777 if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
1781 else printf("O.K.\n");
1783 /*=============================================*/
1785 /*=============================================*/
1789 static char *msg[] = {
1790 "FAILUREs encountered =",
1791 "SERIOUS DEFECTs discovered =",
1792 "DEFECTs discovered =",
1793 "FLAWs discovered =" };
1795 for(i = 0; i < 4; i++) if (ErrCnt[i])
1796 printf("The number of %-29s %d.\n",
1800 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
1801 + ErrCnt[Flaw]) > 0) {
1802 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
1803 Defect] == 0) && (ErrCnt[Flaw] > 0)) {
1804 printf("The arithmetic diagnosed seems ");
1805 printf("Satisfactory though flawed.\n");
1807 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
1808 && ( ErrCnt[Defect] > 0)) {
1809 printf("The arithmetic diagnosed may be Acceptable\n");
1810 printf("despite inconvenient Defects.\n");
1812 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
1813 printf("The arithmetic diagnosed has ");
1814 printf("unacceptable Serious Defects.\n");
1816 if (ErrCnt[Failure] > 0) {
1817 printf("Potentially fatal FAILURE may have spoiled this");
1818 printf(" program's subsequent diagnoses.\n");
1822 printf("No failures, defects nor flaws have been discovered.\n");
1823 if (! ((RMult == Rounded) && (RDiv == Rounded)
1824 && (RAddSub == Rounded) && (RSqrt == Rounded)))
1825 printf("The arithmetic diagnosed seems Satisfactory.\n");
1827 if (StickyBit >= One &&
1828 (Radix - Two) * (Radix - Nine - One) == Zero) {
1829 printf("Rounding appears to conform to ");
1830 printf("the proposed IEEE standard P");
1831 if ((Radix == Two) &&
1832 ((Precision - Four * Three * Two) *
1833 ( Precision - TwentySeven -
1834 TwentySeven + One) == Zero))
1837 if (IEEE) printf(".\n");
1839 printf(",\nexcept for possibly Double Rounding");
1840 printf(" during Gradual Underflow.\n");
1843 printf("The arithmetic diagnosed appears to be Excellent!\n");
1847 printf("\nA total of %d floating point exceptions were registered.\n",
1849 printf("END OF TEST.\n");
1854 #include "paranoia.h"
1861 { return X >= 0. ? 1.0 : -1.0; }
1870 printf("\nTo continue, press RETURN");
1872 read(KEYBOARD, ch, 8);
1874 printf("\nDiagnosis resumes after milestone Number %d", Milestone);
1875 printf(" Page: %d\n\n", PageNo);
1882 TstCond (K, Valid, T)
1885 { if (! Valid) { BadCond(K,T); printf(".\n"); } }
1891 static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
1893 ErrCnt [K] = ErrCnt [K] + 1;
1894 printf("%s: %s", msg[K], T);
1899 X = (Random1 + Random9)^5
1900 Random1 = X - FLOOR(X) + 0.000005 * X;
1901 and returns the new value of Random1
1908 X = Random1 + Random9;
1913 Random1 = Y + X * 0.000005;
1926 SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
1928 if (SqEr < MinSqEr) MinSqEr = SqEr;
1929 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1931 BadCond(ErrKind, "\n");
1932 printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
1933 printf("\tinstead of correct value 0 .\n");
1942 X = FLOOR(Half - X / Radix) * Radix + X;
1943 Q = (Q - X * Z) / Radix + X * X * (D / Radix);
1944 Z = Z - Two * X * D;
1956 if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
1959 Y2 = (X2 - Z2) - (Y - Z2);
1960 X2 = X8 / (Y - Half);
1961 X2 = X2 - Half * X2 * X2;
1962 SqEr = (Y2 + Half) + (Half - X2);
1963 if (SqEr < MinSqEr) MinSqEr = SqEr;
1965 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1975 if (Z == Zero && Q <= Zero)
1976 printf("WARNING: computing\n");
1977 else BadCond(Defect, "computing\n");
1978 printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
1979 printf("\tyielded %.17e;\n", Y);
1980 printf("\twhich compared unequal to correct %.17e ;\n",
1982 printf("\t\tthey differ by %.17e .\n", Y - X);
1984 N = N + 1; /* ... count discrepancies. */
2001 /* PrintIfNPositive */
2005 if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
2014 printf("Since comparison denies Z = 0, evaluating ");
2015 printf("(Z + Z) / Z should be safe.\n");
2017 if (setjmp(ovfl_buf)) goto very_serious;
2019 printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
2021 if (FABS(Q9 - Two) < Radix * U2) {
2022 printf("This is O.K., provided Over/Underflow");
2023 printf(" has NOT just been signaled.\n");
2026 if ((Q9 < One) || (Q9 > Two)) {
2029 ErrCnt [Serious] = ErrCnt [Serious] + 1;
2030 printf("This is a VERY SERIOUS DEFECT!\n");
2034 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2035 printf("This is a DEFECT!\n");
2044 if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
2049 BadCond(Defect, "What prints as Z = ");
2050 printf("%.17e\n\tcompares different from ", Z);
2051 if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
2052 if (! ((Z == Random2)
2053 || (Random2 == Random1)))
2054 printf("1 * Z == %g\n", Random2);
2055 if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
2056 if (Random2 != Random1) {
2057 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2058 BadCond(Defect, "Multiplication does not commute!\n");
2059 printf("\tComparison alleges that 1 * Z = %.17e\n",
2061 printf("\tdiffers from Z * 1 = %.17e\n", Random1);
2071 printf("%s test appears to be inconsistent...\n", s);
2072 printf(" PLEASE NOTIFY KARPINKSI!\n");
2081 { while(*s) printf("%s\n", *s++); }
2085 static char *instr[] = {
2086 "Lest this program stop prematurely, i.e. before displaying\n",
2087 " `END OF TEST',\n",
2088 "try to persuade the computer NOT to terminate execution when an",
2089 "error like Over/Underflow or Division by Zero occurs, but rather",
2090 "to persevere with a surrogate value after, perhaps, displaying some",
2091 "warning. If persuasion avails naught, don't despair but run this",
2092 "program anyway to see how many milestones it passes, and then",
2093 "amend it to make further progress.\n",
2094 "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
2104 static char *head[] = {
2105 "Users are invited to help debug and augment this program so it will",
2106 "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
2107 "Please send suggestions and interesting results to",
2108 "\tRichard Karpinski",
2109 "\tComputer Center U-76",
2110 "\tUniversity of California",
2111 "\tSan Francisco, CA 94143-0704, USA\n",
2112 "In doing so, please include the following information:",
2114 "\tPrecision:\tsingle;",
2116 "\tPrecision:\tdouble;",
2118 "\tVersion:\t10 February 1989;",
2121 "\tOptimization level:\n",
2122 "\tOther relevant compiler options:",
2128 /* Characteristics */
2132 static char *chars[] = {
2133 "Running this program should reveal these characteristics:",
2134 " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
2135 " Precision = number of significant digits carried.",
2136 " U2 = Radix/Radix^Precision = One Ulp",
2137 "\t(OneUlpnit in the Last Place) of 1.000xxx .",
2138 " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
2139 " Adequacy of guard digits for Mult., Div. and Subt.",
2140 " Whether arithmetic is chopped, correctly rounded, or something else",
2141 "\tfor Mult., Div., Add/Subt. and Sqrt.",
2142 " Whether a Sticky Bit used correctly for rounding.",
2143 " UnderflowThreshold = an underflow threshold.",
2144 " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
2145 " V = an overflow threshold, roughly.",
2146 " V0 tells, roughly, whether Infinity is represented.",
2147 " Comparisions are checked for consistency with subtraction",
2148 "\tand for contamination with pseudo-zeros.",
2149 " Sqrt is tested. Y^X is not tested.",
2150 " Extra-precise subexpressions are revealed but NOT YET tested.",
2151 " Decimal-Binary conversion is NOT YET tested for accuracy.",
2160 /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
2161 with further massaging by David M. Gay. */
2163 static char *hist[] = {
2164 "The program attempts to discriminate among",
2165 " FLAWs, like lack of a sticky bit,",
2166 " Serious DEFECTs, like lack of a guard digit, and",
2167 " FAILUREs, like 2+2 == 5 .",
2168 "Failures may confound subsequent diagnoses.\n",
2169 "The diagnostic capabilities of this program go beyond an earlier",
2170 "program called `MACHAR', which can be found at the end of the",
2171 "book `Software Manual for the Elementary Functions' (1980) by",
2172 "W. J. Cody and W. Waite. Although both programs try to discover",
2173 "the Radix, Precision and range (over/underflow thresholds)",
2174 "of the arithmetic, this program tries to cope with a wider variety",
2175 "of pathologies, and to say how well the arithmetic is implemented.",
2176 "\nThe program is based upon a conventional radix representation for",
2177 "floating-point numbers, but also allows logarithmic encoding",
2178 "as used by certain early WANG machines.\n",
2179 "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
2180 "see source comments for more history.",
2187 pow(x, y) /* return x ^ y (exponentiation) */
2190 extern double exp(), frexp(), ldexp(), log(), modf();
2193 int ex, ey = 0, flip = 0;
2197 if ((y < -1100. || y > 1100.) && x != -1.) return exp(y * log(x));
2199 if (y < 0.) { y = -y; flip = 1; }
2201 if (y) xy = exp(y * log(x));
2203 /* next several lines assume >= 32 bit integers */
2205 if (i = ye) for(;;) {
2206 if (i & 1) { xy *= x; ey += ex; }
2207 if (!(i >>= 1)) break;
2210 if (x < .5) { x *= 2.; ex -= 1; }
2212 if (flip) { xy = 1. / xy; ey = -ey; }
2213 return ldexp(xy, ey);
2216 #endif /* NO_FLOATS */