diff options
author | zakk <zakk@edf5b092-35ff-0310-97b2-ce42778d08ea> | 2005-08-26 17:39:27 +0000 |
---|---|---|
committer | zakk <zakk@edf5b092-35ff-0310-97b2-ce42778d08ea> | 2005-08-26 17:39:27 +0000 |
commit | 6bf20c78f5b69d40bcc4931df93d29198435ab67 (patch) | |
tree | e3eda937a05d7db42de725b7013bd0344b987f34 /lcc/tst/paranoia.c | |
parent | 872d4d7f55af706737ffb361bb76ad13e7496770 (diff) | |
download | ioquake3-aero-6bf20c78f5b69d40bcc4931df93d29198435ab67.tar.gz ioquake3-aero-6bf20c78f5b69d40bcc4931df93d29198435ab67.zip |
newlines fixed
git-svn-id: svn://svn.icculus.org/quake3/trunk@6 edf5b092-35ff-0310-97b2-ce42778d08ea
Diffstat (limited to 'lcc/tst/paranoia.c')
-rwxr-xr-x | lcc/tst/paranoia.c | 4406 |
1 files changed, 2203 insertions, 2203 deletions
diff --git a/lcc/tst/paranoia.c b/lcc/tst/paranoia.c index 76e0f05..ec9f8ce 100755 --- a/lcc/tst/paranoia.c +++ b/lcc/tst/paranoia.c @@ -1,2203 +1,2203 @@ -#undef V9
-#define NOPAUSE
-/* A C version of Kahan's Floating Point Test "Paranoia"
-
- Thos Sumner, UCSF, Feb. 1985
- David Gay, BTL, Jan. 1986
-
- This is a rewrite from the Pascal version by
-
- B. A. Wichmann, 18 Jan. 1985
-
- (and does NOT exhibit good C programming style).
-
-(C) Apr 19 1983 in BASIC version by:
- Professor W. M. Kahan,
- 567 Evans Hall
- Electrical Engineering & Computer Science Dept.
- University of California
- Berkeley, California 94720
- USA
-
-converted to Pascal by:
- B. A. Wichmann
- National Physical Laboratory
- Teddington Middx
- TW11 OLW
- UK
-
-converted to C by:
-
- David M. Gay and Thos Sumner
- AT&T Bell Labs Computer Center, Rm. U-76
- 600 Mountain Avenue University of California
- Murray Hill, NJ 07974 San Francisco, CA 94143
- USA USA
-
-with simultaneous corrections to the Pascal source (reflected
-in the Pascal source available over netlib).
-[A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
-
-Reports of results on various systems from all the versions
-of Paranoia are being collected by Richard Karpinski at the
-same address as Thos Sumner. This includes sample outputs,
-bug reports, and criticisms.
-
-You may copy this program freely if you acknowledge its source.
-Comments on the Pascal version to NPL, please.
-
-
-The C version catches signals from floating-point exceptions.
-If signal(SIGFPE,...) is unavailable in your environment, you may
-#define NOSIGNAL to comment out the invocations of signal.
-
-This source file is too big for some C compilers, but may be split
-into pieces. Comments containing "SPLIT" suggest convenient places
-for this splitting. At the end of these comments is an "ed script"
-(for the UNIX(tm) editor ed) that will do this splitting.
-
-By #defining Single when you compile this source, you may obtain
-a single-precision C version of Paranoia.
-
-
-The following is from the introductory commentary from Wichmann's work:
-
-The BASIC program of Kahan is written in Microsoft BASIC using many
-facilities which have no exact analogy in Pascal. The Pascal
-version below cannot therefore be exactly the same. Rather than be
-a minimal transcription of the BASIC program, the Pascal coding
-follows the conventional style of block-structured languages. Hence
-the Pascal version could be useful in producing versions in other
-structured languages.
-
-Rather than use identifiers of minimal length (which therefore have
-little mnemonic significance), the Pascal version uses meaningful
-identifiers as follows [Note: A few changes have been made for C]:
-
-
-BASIC C BASIC C BASIC C
-
- A J S StickyBit
- A1 AInverse J0 NoErrors T
- B Radix [Failure] T0 Underflow
- B1 BInverse J1 NoErrors T2 ThirtyTwo
- B2 RadixD2 [SeriousDefect] T5 OneAndHalf
- B9 BMinusU2 J2 NoErrors T7 TwentySeven
- C [Defect] T8 TwoForty
- C1 CInverse J3 NoErrors U OneUlp
- D [Flaw] U0 UnderflowThreshold
- D4 FourD K PageNo U1
- E0 L Milestone U2
- E1 M V
- E2 Exp2 N V0
- E3 N1 V8
- E5 MinSqEr O Zero V9
- E6 SqEr O1 One W
- E7 MaxSqEr O2 Two X
- E8 O3 Three X1
- E9 O4 Four X8
- F1 MinusOne O5 Five X9 Random1
- F2 Half O8 Eight Y
- F3 Third O9 Nine Y1
- F6 P Precision Y2
- F9 Q Y9 Random2
- G1 GMult Q8 Z
- G2 GDiv Q9 Z0 PseudoZero
- G3 GAddSub R Z1
- H R1 RMult Z2
- H1 HInverse R2 RDiv Z9
- I R3 RAddSub
- IO NoTrials R4 RSqrt
- I3 IEEE R9 Random9
-
- SqRWrng
-
-All the variables in BASIC are true variables and in consequence,
-the program is more difficult to follow since the "constants" must
-be determined (the glossary is very helpful). The Pascal version
-uses Real constants, but checks are added to ensure that the values
-are correctly converted by the compiler.
-
-The major textual change to the Pascal version apart from the
-identifiersis that named procedures are used, inserting parameters
-wherehelpful. New procedures are also introduced. The
-correspondence is as follows:
-
-
-BASIC Pascal
-lines
-
- 90- 140 Pause
- 170- 250 Instructions
- 380- 460 Heading
- 480- 670 Characteristics
- 690- 870 History
-2940-2950 Random
-3710-3740 NewD
-4040-4080 DoesYequalX
-4090-4110 PrintIfNPositive
-4640-4850 TestPartialUnderflow
-
-=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
-
-Below is an "ed script" that splits para.c into 10 files
-of the form part[1-8].c, subs.c, and msgs.c, plus a header
-file, paranoia.h, that these files require.
-
-r paranoia.c
-$
-?SPLIT
-+,$w msgs.c
- .,$d
-?SPLIT
- .d
-+d
--,$w subs.c
--,$d
-?part8
-+d
-?include
- .,$w part8.c
- .,$d
--d
-?part7
-+d
-?include
- .,$w part7.c
- .,$d
--d
-?part6
-+d
-?include
- .,$w part6.c
- .,$d
--d
-?part5
-+d
-?include
- .,$w part5.c
- .,$d
--d
-?part4
-+d
-?include
- .,$w part4.c
- .,$d
--d
-?part3
-+d
-?include
- .,$w part3.c
- .,$d
--d
-?part2
-+d
-?include
- .,$w part2.c
- .,$d
-?SPLIT
- .d
-1,/^#include/-1d
-1,$w part1.c
-/Computed constants/,$d
-1,$s/^int/extern &/
-1,$s/^FLOAT/extern &/
-1,$s/^char/extern &/
-1,$s! = .*!;!
-/^Guard/,/^Round/s/^/extern /
-/^jmp_buf/s/^/extern /
-/^Sig_type/s/^/extern /
-s/$/\
-extern void sigfpe();/
-w paranoia.h
-q
-
-*/
-
-#include <stdio.h>
-#ifndef NOSIGNAL
-#include <signal.h>
-#endif
-#include <setjmp.h>
-
-extern double fabs(), floor(), log(), pow(), sqrt();
-
-#ifdef Single
-#define FLOAT float
-#define FABS(x) (float)fabs((double)(x))
-#define FLOOR(x) (float)floor((double)(x))
-#define LOG(x) (float)log((double)(x))
-#define POW(x,y) (float)pow((double)(x),(double)(y))
-#define SQRT(x) (float)sqrt((double)(x))
-#else
-#define FLOAT double
-#define FABS(x) fabs(x)
-#define FLOOR(x) floor(x)
-#define LOG(x) log(x)
-#define POW(x,y) pow(x,y)
-#define SQRT(x) sqrt(x)
-#endif
-
-jmp_buf ovfl_buf;
-typedef void (*Sig_type)();
-Sig_type sigsave;
-
-#define KEYBOARD 0
-
-FLOAT Radix, BInvrse, RadixD2, BMinusU2;
-FLOAT Sign(), Random();
-
-/*Small floating point constants.*/
-FLOAT Zero = 0.0;
-FLOAT Half = 0.5;
-FLOAT One = 1.0;
-FLOAT Two = 2.0;
-FLOAT Three = 3.0;
-FLOAT Four = 4.0;
-FLOAT Five = 5.0;
-FLOAT Eight = 8.0;
-FLOAT Nine = 9.0;
-FLOAT TwentySeven = 27.0;
-FLOAT ThirtyTwo = 32.0;
-FLOAT TwoForty = 240.0;
-FLOAT MinusOne = -1.0;
-FLOAT OneAndHalf = 1.5;
-/*Integer constants*/
-int NoTrials = 20; /*Number of tests for commutativity. */
-#define False 0
-#define True 1
-
-/* Definitions for declared types
- Guard == (Yes, No);
- Rounding == (Chopped, Rounded, Other);
- Message == packed array [1..40] of char;
- Class == (Flaw, Defect, Serious, Failure);
- */
-#define Yes 1
-#define No 0
-#define Chopped 2
-#define Rounded 1
-#define Other 0
-#define Flaw 3
-#define Defect 2
-#define Serious 1
-#define Failure 0
-typedef int Guard, Rounding, Class;
-typedef char Message;
-
-/* Declarations of Variables */
-int Indx;
-char ch[8];
-FLOAT AInvrse, A1;
-FLOAT C, CInvrse;
-FLOAT D, FourD;
-FLOAT E0, E1, Exp2, E3, MinSqEr;
-FLOAT SqEr, MaxSqEr, E9;
-FLOAT Third;
-FLOAT F6, F9;
-FLOAT H, HInvrse;
-int I;
-FLOAT StickyBit, J;
-FLOAT MyZero;
-FLOAT Precision;
-FLOAT Q, Q9;
-FLOAT R, Random9;
-FLOAT T, Underflow, S;
-FLOAT OneUlp, UfThold, U1, U2;
-FLOAT V, V0, V9;
-FLOAT W;
-FLOAT X, X1, X2, X8, Random1;
-FLOAT Y, Y1, Y2, Random2;
-FLOAT Z, PseudoZero, Z1, Z2, Z9;
-int ErrCnt[4];
-int fpecount;
-int Milestone;
-int PageNo;
-int M, N, N1;
-Guard GMult, GDiv, GAddSub;
-Rounding RMult, RDiv, RAddSub, RSqrt;
-int Break, Done, NotMonot, Monot, Anomaly, IEEE,
- SqRWrng, UfNGrad;
-/* Computed constants. */
-/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
-/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
-
-/* floating point exception receiver */
- void
-sigfpe(i)
-{
- fpecount++;
- printf("\n* * * FLOATING-POINT ERROR * * *\n");
- fflush(stdout);
- if (sigsave) {
-#ifndef NOSIGNAL
- signal(SIGFPE, sigsave);
-#endif
- sigsave = 0;
- longjmp(ovfl_buf, 1);
- }
- abort();
-}
-
-main()
-{
-#ifdef mc
- char *out;
- ieee_flags("set", "precision", "double", &out);
-#endif
- /* First two assignments use integer right-hand sides. */
- Zero = 0;
- One = 1;
- Two = One + One;
- Three = Two + One;
- Four = Three + One;
- Five = Four + One;
- Eight = Four + Four;
- Nine = Three * Three;
- TwentySeven = Nine * Three;
- ThirtyTwo = Four * Eight;
- TwoForty = Four * Five * Three * Four;
- MinusOne = -One;
- Half = One / Two;
- OneAndHalf = One + Half;
- ErrCnt[Failure] = 0;
- ErrCnt[Serious] = 0;
- ErrCnt[Defect] = 0;
- ErrCnt[Flaw] = 0;
- PageNo = 1;
- /*=============================================*/
- Milestone = 0;
- /*=============================================*/
-#ifndef NOSIGNAL
- signal(SIGFPE, sigfpe);
-#endif
- Instructions();
- Pause();
- Heading();
- Pause();
- Characteristics();
- Pause();
- History();
- Pause();
- /*=============================================*/
- Milestone = 7;
- /*=============================================*/
- printf("Program is now RUNNING tests on small integers:\n");
-
- TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
- && (One > Zero) && (One + One == Two),
- "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
- Z = - Zero;
- if (Z != 0.0) {
- ErrCnt[Failure] = ErrCnt[Failure] + 1;
- printf("Comparison alleges that -0.0 is Non-zero!\n");
- U1 = 0.001;
- Radix = 1;
- TstPtUf();
- }
- TstCond (Failure, (Three == Two + One) && (Four == Three + One)
- && (Four + Two * (- Two) == Zero)
- && (Four - Three - One == Zero),
- "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
- TstCond (Failure, (MinusOne == (0 - One))
- && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
- && (MinusOne + FABS(One) == Zero)
- && (MinusOne + MinusOne * MinusOne == Zero),
- "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
- TstCond (Failure, Half + MinusOne + Half == Zero,
- "1/2 + (-1) + 1/2 != 0");
- /*=============================================*/
- /*SPLIT
- part2();
- part3();
- part4();
- part5();
- part6();
- part7();
- part8();
- }
-#include "paranoia.h"
-part2(){
-*/
- Milestone = 10;
- /*=============================================*/
- TstCond (Failure, (Nine == Three * Three)
- && (TwentySeven == Nine * Three) && (Eight == Four + Four)
- && (ThirtyTwo == Eight * Four)
- && (ThirtyTwo - TwentySeven - Four - One == Zero),
- "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
- TstCond (Failure, (Five == Four + One) &&
- (TwoForty == Four * Five * Three * Four)
- && (TwoForty / Three - Four * Four * Five == Zero)
- && ( TwoForty / Four - Five * Three * Four == Zero)
- && ( TwoForty / Five - Four * Three * Four == Zero),
- "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
- if (ErrCnt[Failure] == 0) {
- printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
- printf("\n");
- }
- printf("Searching for Radix and Precision.\n");
- W = One;
- do {
- W = W + W;
- Y = W + One;
- Z = Y - W;
- Y = Z - One;
- } while (MinusOne + FABS(Y) < Zero);
- /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
- Precision = Zero;
- Y = One;
- do {
- Radix = W + Y;
- Y = Y + Y;
- Radix = Radix - W;
- } while ( Radix == Zero);
- if (Radix < Two) Radix = One;
- printf("Radix = %f .\n", Radix);
- if (Radix != 1) {
- W = One;
- do {
- Precision = Precision + One;
- W = W * Radix;
- Y = W + One;
- } while ((Y - W) == One);
- }
- /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
- ...*/
- U1 = One / W;
- U2 = Radix * U1;
- printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
- printf("Recalculating radix and precision\n ");
-
- /*save old values*/
- E0 = Radix;
- E1 = U1;
- E9 = U2;
- E3 = Precision;
-
- X = Four / Three;
- Third = X - One;
- F6 = Half - Third;
- X = F6 + F6;
- X = FABS(X - Third);
- if (X < U2) X = U2;
-
- /*... now X = (unknown no.) ulps of 1+...*/
- do {
- U2 = X;
- Y = Half * U2 + ThirtyTwo * U2 * U2;
- Y = One + Y;
- X = Y - One;
- } while ( ! ((U2 <= X) || (X <= Zero)));
-
- /*... now U2 == 1 ulp of 1 + ... */
- X = Two / Three;
- F6 = X - Half;
- Third = F6 + F6;
- X = Third - Half;
- X = FABS(X + F6);
- if (X < U1) X = U1;
-
- /*... now X == (unknown no.) ulps of 1 -... */
- do {
- U1 = X;
- Y = Half * U1 + ThirtyTwo * U1 * U1;
- Y = Half - Y;
- X = Half + Y;
- Y = Half - X;
- X = Half + Y;
- } while ( ! ((U1 <= X) || (X <= Zero)));
- /*... now U1 == 1 ulp of 1 - ... */
- if (U1 == E1) printf("confirms closest relative separation U1 .\n");
- else printf("gets better closest relative separation U1 = %.7e .\n", U1);
- W = One / U1;
- F9 = (Half - U1) + Half;
- Radix = FLOOR(0.01 + U2 / U1);
- if (Radix == E0) printf("Radix confirmed.\n");
- else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
- TstCond (Defect, Radix <= Eight + Eight,
- "Radix is too big: roundoff problems");
- TstCond (Flaw, (Radix == Two) || (Radix == 10)
- || (Radix == One), "Radix is not as good as 2 or 10");
- /*=============================================*/
- Milestone = 20;
- /*=============================================*/
- TstCond (Failure, F9 - Half < Half,
- "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
- X = F9;
- I = 1;
- Y = X - Half;
- Z = Y - Half;
- TstCond (Failure, (X != One)
- || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
- X = One + U2;
- I = 0;
- /*=============================================*/
- Milestone = 25;
- /*=============================================*/
- /*... BMinusU2 = nextafter(Radix, 0) */
- BMinusU2 = Radix - One;
- BMinusU2 = (BMinusU2 - U2) + One;
- /* Purify Integers */
- if (Radix != One) {
- X = - TwoForty * LOG(U1) / LOG(Radix);
- Y = FLOOR(Half + X);
- if (FABS(X - Y) * Four < One) X = Y;
- Precision = X / TwoForty;
- Y = FLOOR(Half + Precision);
- if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
- }
- if ((Precision != FLOOR(Precision)) || (Radix == One)) {
- printf("Precision cannot be characterized by an Integer number\n");
- printf("of significant digits but, by itself, this is a minor flaw.\n");
- }
- if (Radix == One)
- printf("logarithmic encoding has precision characterized solely by U1.\n");
- else printf("The number of significant digits of the Radix is %f .\n",
- Precision);
- TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
- "Precision worse than 5 decimal figures ");
- /*=============================================*/
- Milestone = 30;
- /*=============================================*/
- /* Test for extra-precise subepressions */
- X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
- do {
- Z2 = X;
- X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
- } while ( ! ((Z2 <= X) || (X <= Zero)));
- X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
- do {
- Z1 = Z;
- Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
- + One / Two)) + One / Two;
- } while ( ! ((Z1 <= Z) || (Z <= Zero)));
- do {
- do {
- Y1 = Y;
- Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
- )) + Half;
- } while ( ! ((Y1 <= Y) || (Y <= Zero)));
- X1 = X;
- X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
- } while ( ! ((X1 <= X) || (X <= Zero)));
- if ((X1 != Y1) || (X1 != Z1)) {
- BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
- printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1);
- printf("are symptoms of inconsistencies introduced\n");
- printf("by extra-precise evaluation of arithmetic subexpressions.\n");
- notify("Possibly some part of this");
- if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf(
- "That feature is not tested further by this program.\n") ;
- }
- else {
- if ((Z1 != U1) || (Z2 != U2)) {
- if ((Z1 >= U1) || (Z2 >= U2)) {
- BadCond(Failure, "");
- notify("Precision");
- printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
- printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
- }
- else {
- if ((Z1 <= Zero) || (Z2 <= Zero)) {
- printf("Because of unusual Radix = %f", Radix);
- printf(", or exact rational arithmetic a result\n");
- printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
- notify("of an\nextra-precision");
- }
- if (Z1 != Z2 || Z1 > Zero) {
- X = Z1 / U1;
- Y = Z2 / U2;
- if (Y > X) X = Y;
- Q = - LOG(X);
- printf("Some subexpressions appear to be calculated extra\n");
- printf("precisely with about %g extra B-digits, i.e.\n",
- (Q / LOG(Radix)));
- printf("roughly %g extra significant decimals.\n",
- Q / LOG(10.));
- }
- printf("That feature is not tested further by this program.\n");
- }
- }
- }
- Pause();
- /*=============================================*/
- /*SPLIT
- }
-#include "paranoia.h"
-part3(){
-*/
- Milestone = 35;
- /*=============================================*/
- if (Radix >= Two) {
- X = W / (Radix * Radix);
- Y = X + One;
- Z = Y - X;
- T = Z + U2;
- X = T - Z;
- TstCond (Failure, X == U2,
- "Subtraction is not normalized X=Y,X+Z != Y+Z!");
- if (X == U2) printf(
- "Subtraction appears to be normalized, as it should be.");
- }
- printf("\nChecking for guard digit in *, /, and -.\n");
- Y = F9 * One;
- Z = One * F9;
- X = F9 - Half;
- Y = (Y - Half) - X;
- Z = (Z - Half) - X;
- X = One + U2;
- T = X * Radix;
- R = Radix * X;
- X = T - Radix;
- X = X - Radix * U2;
- T = R - Radix;
- T = T - Radix * U2;
- X = X * (Radix - One);
- T = T * (Radix - One);
- if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
- else {
- GMult = No;
- TstCond (Serious, False,
- "* lacks a Guard Digit, so 1*X != X");
- }
- Z = Radix * U2;
- X = One + Z;
- Y = FABS((X + Z) - X * X) - U2;
- X = One - U2;
- Z = FABS((X - U2) - X * X) - U1;
- TstCond (Failure, (Y <= Zero)
- && (Z <= Zero), "* gets too many final digits wrong.\n");
- Y = One - U2;
- X = One + U2;
- Z = One / Y;
- Y = Z - X;
- X = One / Three;
- Z = Three / Nine;
- X = X - Z;
- T = Nine / TwentySeven;
- Z = Z - T;
- TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
- "Division lacks a Guard Digit, so error can exceed 1 ulp\nor 1/3 and 3/9 and 9/27 may disagree");
- Y = F9 / One;
- X = F9 - Half;
- Y = (Y - Half) - X;
- X = One + U2;
- T = X / One;
- X = T - X;
- if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
- else {
- GDiv = No;
- TstCond (Serious, False,
- "Division lacks a Guard Digit, so X/1 != X");
- }
- X = One / (One + U2);
- Y = X - Half - Half;
- TstCond (Serious, Y < Zero,
- "Computed value of 1/1.000..1 >= 1");
- X = One - U2;
- Y = One + Radix * U2;
- Z = X * Radix;
- T = Y * Radix;
- R = Z / Radix;
- StickyBit = T / Radix;
- X = R - X;
- Y = StickyBit - Y;
- TstCond (Failure, X == Zero && Y == Zero,
- "* and/or / gets too many last digits wrong");
- Y = One - U1;
- X = One - F9;
- Y = One - Y;
- T = Radix - U2;
- Z = Radix - BMinusU2;
- T = Radix - T;
- if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
- else {
- GAddSub = No;
- TstCond (Serious, False,
- "- lacks Guard Digit, so cancellation is obscured");
- }
- if (F9 != One && F9 - One >= Zero) {
- BadCond(Serious, "comparison alleges (1-U1) < 1 although\n");
- printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n");
- printf(" such precautions against division by zero as\n");
- printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
- }
- if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
- " *, /, and - appear to have guard digits, as they should.\n");
- /*=============================================*/
- Milestone = 40;
- /*=============================================*/
- Pause();
- printf("Checking rounding on multiply, divide and add/subtract.\n");
- RMult = Other;
- RDiv = Other;
- RAddSub = Other;
- RadixD2 = Radix / Two;
- A1 = Two;
- Done = False;
- do {
- AInvrse = Radix;
- do {
- X = AInvrse;
- AInvrse = AInvrse / A1;
- } while ( ! (FLOOR(AInvrse) != AInvrse));
- Done = (X == One) || (A1 > Three);
- if (! Done) A1 = Nine + One;
- } while ( ! (Done));
- if (X == One) A1 = Radix;
- AInvrse = One / A1;
- X = A1;
- Y = AInvrse;
- Done = False;
- do {
- Z = X * Y - Half;
- TstCond (Failure, Z == Half,
- "X * (1/X) differs from 1");
- Done = X == Radix;
- X = Radix;
- Y = One / X;
- } while ( ! (Done));
- Y2 = One + U2;
- Y1 = One - U2;
- X = OneAndHalf - U2;
- Y = OneAndHalf + U2;
- Z = (X - U2) * Y2;
- T = Y * Y1;
- Z = Z - X;
- T = T - X;
- X = X * Y2;
- Y = (Y + U2) * Y1;
- X = X - OneAndHalf;
- Y = Y - OneAndHalf;
- if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
- X = (OneAndHalf + U2) * Y2;
- Y = OneAndHalf - U2 - U2;
- Z = OneAndHalf + U2 + U2;
- T = (OneAndHalf - U2) * Y1;
- X = X - (Z + U2);
- StickyBit = Y * Y1;
- S = Z * Y2;
- T = T - Y;
- Y = (U2 - Y) + StickyBit;
- Z = S - (Z + U2 + U2);
- StickyBit = (Y2 + U2) * Y1;
- Y1 = Y2 * Y1;
- StickyBit = StickyBit - Y2;
- Y1 = Y1 - Half;
- if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
- && ( StickyBit == Zero) && (Y1 == Half)) {
- RMult = Rounded;
- printf("Multiplication appears to round correctly.\n");
- }
- else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
- && (T < Zero) && (StickyBit + U2 == Zero)
- && (Y1 < Half)) {
- RMult = Chopped;
- printf("Multiplication appears to chop.\n");
- }
- else printf("* is neither chopped nor correctly rounded.\n");
- if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
- }
- else printf("* is neither chopped nor correctly rounded.\n");
- /*=============================================*/
- Milestone = 45;
- /*=============================================*/
- Y2 = One + U2;
- Y1 = One - U2;
- Z = OneAndHalf + U2 + U2;
- X = Z / Y2;
- T = OneAndHalf - U2 - U2;
- Y = (T - U2) / Y1;
- Z = (Z + U2) / Y2;
- X = X - OneAndHalf;
- Y = Y - T;
- T = T / Y1;
- Z = Z - (OneAndHalf + U2);
- T = (U2 - OneAndHalf) + T;
- if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
- X = OneAndHalf / Y2;
- Y = OneAndHalf - U2;
- Z = OneAndHalf + U2;
- X = X - Y;
- T = OneAndHalf / Y1;
- Y = Y / Y1;
- T = T - (Z + U2);
- Y = Y - Z;
- Z = Z / Y2;
- Y1 = (Y2 + U2) / Y2;
- Z = Z - OneAndHalf;
- Y2 = Y1 - Y2;
- Y1 = (F9 - U1) / F9;
- if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
- && (Y2 == Zero) && (Y2 == Zero)
- && (Y1 - Half == F9 - Half )) {
- RDiv = Rounded;
- printf("Division appears to round correctly.\n");
- if (GDiv == No) notify("Division");
- }
- else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
- && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
- RDiv = Chopped;
- printf("Division appears to chop.\n");
- }
- }
- if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
- BInvrse = One / Radix;
- TstCond (Failure, (BInvrse * Radix - Half == Half),
- "Radix * ( 1 / Radix ) differs from 1");
- /*=============================================*/
- /*SPLIT
- }
-#include "paranoia.h"
-part4(){
-*/
- Milestone = 50;
- /*=============================================*/
- TstCond (Failure, ((F9 + U1) - Half == Half)
- && ((BMinusU2 + U2 ) - One == Radix - One),
- "Incomplete carry-propagation in Addition");
- X = One - U1 * U1;
- Y = One + U2 * (One - U2);
- Z = F9 - Half;
- X = (X - Half) - Z;
- Y = Y - One;
- if ((X == Zero) && (Y == Zero)) {
- RAddSub = Chopped;
- printf("Add/Subtract appears to be chopped.\n");
- }
- if (GAddSub == Yes) {
- X = (Half + U2) * U2;
- Y = (Half - U2) * U2;
- X = One + X;
- Y = One + Y;
- X = (One + U2) - X;
- Y = One - Y;
- if ((X == Zero) && (Y == Zero)) {
- X = (Half + U2) * U1;
- Y = (Half - U2) * U1;
- X = One - X;
- Y = One - Y;
- X = F9 - X;
- Y = One - Y;
- if ((X == Zero) && (Y == Zero)) {
- RAddSub = Rounded;
- printf("Addition/Subtraction appears to round correctly.\n");
- if (GAddSub == No) notify("Add/Subtract");
- }
- else printf("Addition/Subtraction neither rounds nor chops.\n");
- }
- else printf("Addition/Subtraction neither rounds nor chops.\n");
- }
- else printf("Addition/Subtraction neither rounds nor chops.\n");
- S = One;
- X = One + Half * (One + Half);
- Y = (One + U2) * Half;
- Z = X - Y;
- T = Y - X;
- StickyBit = Z + T;
- if (StickyBit != Zero) {
- S = Zero;
- BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
- }
- StickyBit = Zero;
- if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
- && (RMult == Rounded) && (RDiv == Rounded)
- && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
- printf("Checking for sticky bit.\n");
- X = (Half + U1) * U2;
- Y = Half * U2;
- Z = One + Y;
- T = One + X;
- if ((Z - One <= Zero) && (T - One >= U2)) {
- Z = T + Y;
- Y = Z - X;
- if ((Z - T >= U2) && (Y - T == Zero)) {
- X = (Half + U1) * U1;
- Y = Half * U1;
- Z = One - Y;
- T = One - X;
- if ((Z - One == Zero) && (T - F9 == Zero)) {
- Z = (Half - U1) * U1;
- T = F9 - Z;
- Q = F9 - Y;
- if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
- Z = (One + U2) * OneAndHalf;
- T = (OneAndHalf + U2) - Z + U2;
- X = One + Half / Radix;
- Y = One + Radix * U2;
- Z = X * Y;
- if (T == Zero && X + Radix * U2 - Z == Zero) {
- if (Radix != Two) {
- X = Two + U2;
- Y = X / Two;
- if ((Y - One == Zero)) StickyBit = S;
- }
- else StickyBit = S;
- }
- }
- }
- }
- }
- }
- if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
- else printf("Sticky bit used incorrectly or not at all.\n");
- TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
- RMult == Other || RDiv == Other || RAddSub == Other),
- "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");
- /*=============================================*/
- Milestone = 60;
- /*=============================================*/
- printf("\n");
- printf("Does Multiplication commute? ");
- printf("Testing on %d random pairs.\n", NoTrials);
- Random9 = SQRT(3.0);
- Random1 = Third;
- I = 1;
- do {
- X = Random();
- Y = Random();
- Z9 = Y * X;
- Z = X * Y;
- Z9 = Z - Z9;
- I = I + 1;
- } while ( ! ((I > NoTrials) || (Z9 != Zero)));
- if (I == NoTrials) {
- Random1 = One + Half / Three;
- Random2 = (U2 + U1) + One;
- Z = Random1 * Random2;
- Y = Random2 * Random1;
- Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
- Three) * ((U2 + U1) + One);
- }
- if (! ((I == NoTrials) || (Z9 == Zero)))
- BadCond(Defect, "X * Y == Y * X trial fails.\n");
- else printf(" No failures found in %d integer pairs.\n", NoTrials);
- /*=============================================*/
- Milestone = 70;
- /*=============================================*/
- printf("\nRunning test of square root(x).\n");
- TstCond (Failure, (Zero == SQRT(Zero))
- && (- Zero == SQRT(- Zero))
- && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
- MinSqEr = Zero;
- MaxSqEr = Zero;
- J = Zero;
- X = Radix;
- OneUlp = U2;
- SqXMinX (Serious);
- X = BInvrse;
- OneUlp = BInvrse * U1;
- SqXMinX (Serious);
- X = U1;
- OneUlp = U1 * U1;
- SqXMinX (Serious);
- if (J != Zero) Pause();
- printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
- J = Zero;
- X = Two;
- Y = Radix;
- if ((Radix != One)) do {
- X = Y;
- Y = Radix * Y;
- } while ( ! ((Y - X >= NoTrials)));
- OneUlp = X * U2;
- I = 1;
- while (I <= NoTrials) {
- X = X + One;
- SqXMinX (Defect);
- if (J > Zero) break;
- I = I + 1;
- }
- printf("Test for sqrt monotonicity.\n");
- I = - 1;
- X = BMinusU2;
- Y = Radix;
- Z = Radix + Radix * U2;
- NotMonot = False;
- Monot = False;
- while ( ! (NotMonot || Monot)) {
- I = I + 1;
- X = SQRT(X);
- Q = SQRT(Y);
- Z = SQRT(Z);
- if ((X > Q) || (Q > Z)) NotMonot = True;
- else {
- Q = FLOOR(Q + Half);
- if ((I > 0) || (Radix == Q * Q)) Monot = True;
- else if (I > 0) {
- if (I > 1) Monot = True;
- else {
- Y = Y * BInvrse;
- X = Y - U1;
- Z = Y + U1;
- }
- }
- else {
- Y = Q;
- X = Y - U2;
- Z = Y + U2;
- }
- }
- }
- if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
- else {
- BadCond(Defect, "");
- printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
- }
- /*=============================================*/
- /*SPLIT
- }
-#include "paranoia.h"
-part5(){
-*/
- Milestone = 80;
- /*=============================================*/
- MinSqEr = MinSqEr + Half;
- MaxSqEr = MaxSqEr - Half;
- Y = (SQRT(One + U2) - One) / U2;
- SqEr = (Y - One) + U2 / Eight;
- if (SqEr > MaxSqEr) MaxSqEr = SqEr;
- SqEr = Y + U2 / Eight;
- if (SqEr < MinSqEr) MinSqEr = SqEr;
- Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
- SqEr = Y + U1 / Eight;
- if (SqEr > MaxSqEr) MaxSqEr = SqEr;
- SqEr = (Y + One) + U1 / Eight;
- if (SqEr < MinSqEr) MinSqEr = SqEr;
- OneUlp = U2;
- X = OneUlp;
- for( Indx = 1; Indx <= 3; ++Indx) {
- Y = SQRT((X + U1 + X) + F9);
- Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
- Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
- SqEr = (Y + Half) + Z;
- if (SqEr < MinSqEr) MinSqEr = SqEr;
- SqEr = (Y - Half) + Z;
- if (SqEr > MaxSqEr) MaxSqEr = SqEr;
- if (((Indx == 1) || (Indx == 3)))
- X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
- else {
- OneUlp = U1;
- X = - OneUlp;
- }
- }
- /*=============================================*/
- Milestone = 85;
- /*=============================================*/
- SqRWrng = False;
- Anomaly = False;
- RSqrt = Other; /* ~dgh */
- if (Radix != One) {
- printf("Testing whether sqrt is rounded or chopped.\n");
- D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
- /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
- X = D / Radix;
- Y = D / A1;
- if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
- Anomaly = True;
- }
- else {
- X = Zero;
- Z2 = X;
- Y = One;
- Y2 = Y;
- Z1 = Radix - One;
- FourD = Four * D;
- do {
- if (Y2 > Z2) {
- Q = Radix;
- Y1 = Y;
- do {
- X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
- Q = Y1;
- Y1 = X1;
- } while ( ! (X1 <= Zero));
- if (Q <= One) {
- Z2 = Y2;
- Z = Y;
- }
- }
- Y = Y + Two;
- X = X + Eight;
- Y2 = Y2 + X;
- if (Y2 >= FourD) Y2 = Y2 - FourD;
- } while ( ! (Y >= D));
- X8 = FourD - Z2;
- Q = (X8 + Z * Z) / FourD;
- X8 = X8 / Eight;
- if (Q != FLOOR(Q)) Anomaly = True;
- else {
- Break = False;
- do {
- X = Z1 * Z;
- X = X - FLOOR(X / Radix) * Radix;
- if (X == One)
- Break = True;
- else
- Z1 = Z1 - One;
- } while ( ! (Break || (Z1 <= Zero)));
- if ((Z1 <= Zero) && (! Break)) Anomaly = True;
- else {
- if (Z1 > RadixD2) Z1 = Z1 - Radix;
- do {
- NewD();
- } while ( ! (U2 * D >= F9));
- if (D * Radix - D != W - D) Anomaly = True;
- else {
- Z2 = D;
- I = 0;
- Y = D + (One + Z) * Half;
- X = D + Z + Q;
- SR3750();
- Y = D + (One - Z) * Half + D;
- X = D - Z + D;
- X = X + Q + X;
- SR3750();
- NewD();
- if (D - Z2 != W - Z2) Anomaly = True;
- else {
- Y = (D - Z2) + (Z2 + (One - Z) * Half);
- X = (D - Z2) + (Z2 - Z + Q);
- SR3750();
- Y = (One + Z) * Half;
- X = Q;
- SR3750();
- if (I == 0) Anomaly = True;
- }
- }
- }
- }
- }
- if ((I == 0) || Anomaly) {
- BadCond(Failure, "Anomalous arithmetic with Integer < ");
- printf("Radix^Precision = %.7e\n", W);
- printf(" fails test whether sqrt rounds or chops.\n");
- SqRWrng = True;
- }
- }
- if (! Anomaly) {
- if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
- RSqrt = Rounded;
- printf("Square root appears to be correctly rounded.\n");
- }
- else {
- if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
- || (MinSqEr + Radix < Half)) SqRWrng = True;
- else {
- RSqrt = Chopped;
- printf("Square root appears to be chopped.\n");
- }
- }
- }
- if (SqRWrng) {
- printf("Square root is neither chopped nor correctly rounded.\n");
- printf("Observed errors run from %.7e ", MinSqEr - Half);
- printf("to %.7e ulps.\n", Half + MaxSqEr);
- TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
- "sqrt gets too many last digits wrong");
- }
- /*=============================================*/
- Milestone = 90;
- /*=============================================*/
- Pause();
- printf("Testing powers Z^i for small Integers Z and i.\n");
- N = 0;
- /* ... test powers of zero. */
- I = 0;
- Z = -Zero;
- M = 3.0;
- Break = False;
- do {
- X = One;
- SR3980();
- if (I <= 10) {
- I = 1023;
- SR3980();
- }
- if (Z == MinusOne) Break = True;
- else {
- Z = MinusOne;
- PrintIfNPositive();
- N = 0;
- /* .. if(-1)^N is invalid, replace MinusOne by One. */
- I = - 4;
- }
- } while ( ! Break);
- PrintIfNPositive();
- N1 = N;
- N = 0;
- Z = A1;
- M = FLOOR(Two * LOG(W) / LOG(A1));
- Break = False;
- do {
- X = Z;
- I = 1;
- SR3980();
- if (Z == AInvrse) Break = True;
- else Z = AInvrse;
- } while ( ! (Break));
- /*=============================================*/
- Milestone = 100;
- /*=============================================*/
- /* Powers of Radix have been tested, */
- /* next try a few primes */
- M = NoTrials;
- Z = Three;
- do {
- X = Z;
- I = 1;
- SR3980();
- do {
- Z = Z + Two;
- } while ( Three * FLOOR(Z / Three) == Z );
- } while ( Z < Eight * Three );
- if (N > 0) {
- printf("Errors like this may invalidate financial calculations\n");
- printf("\tinvolving interest rates.\n");
- }
- PrintIfNPositive();
- N += N1;
- if (N == 0) printf("... no discrepancis found.\n");
- if (N > 0) Pause();
- else printf("\n");
- /*=============================================*/
- /*SPLIT
- }
-#include "paranoia.h"
-part6(){
-*/
- Milestone = 110;
- /*=============================================*/
- printf("Seeking Underflow thresholds UfThold and E0.\n");
- D = U1;
- if (Precision != FLOOR(Precision)) {
- D = BInvrse;
- X = Precision;
- do {
- D = D * BInvrse;
- X = X - One;
- } while ( X > Zero);
- }
- Y = One;
- Z = D;
- /* ... D is power of 1/Radix < 1. */
- do {
- C = Y;
- Y = Z;
- Z = Y * Y;
- } while ((Y > Z) && (Z + Z > Z));
- Y = C;
- Z = Y * D;
- do {
- C = Y;
- Y = Z;
- Z = Y * D;
- } while ((Y > Z) && (Z + Z > Z));
- if (Radix < Two) HInvrse = Two;
- else HInvrse = Radix;
- H = One / HInvrse;
- /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
- CInvrse = One / C;
- E0 = C;
- Z = E0 * H;
- /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
- do {
- Y = E0;
- E0 = Z;
- Z = E0 * H;
- } while ((E0 > Z) && (Z + Z > Z));
- UfThold = E0;
- E1 = Zero;
- Q = Zero;
- E9 = U2;
- S = One + E9;
- D = C * S;
- if (D <= C) {
- E9 = Radix * U2;
- S = One + E9;
- D = C * S;
- if (D <= C) {
- BadCond(Failure, "multiplication gets too many last digits wrong.\n");
- Underflow = E0;
- Y1 = Zero;
- PseudoZero = Z;
- Pause();
- }
- }
- else {
- Underflow = D;
- PseudoZero = Underflow * H;
- UfThold = Zero;
- do {
- Y1 = Underflow;
- Underflow = PseudoZero;
- if (E1 + E1 <= E1) {
- Y2 = Underflow * HInvrse;
- E1 = FABS(Y1 - Y2);
- Q = Y1;
- if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
- }
- PseudoZero = PseudoZero * H;
- } while ((Underflow > PseudoZero)
- && (PseudoZero + PseudoZero > PseudoZero));
- }
- /* Comment line 4530 .. 4560 */
- if (PseudoZero != Zero) {
- printf("\n");
- Z = PseudoZero;
- /* ... Test PseudoZero for "phoney- zero" violates */
- /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
- ... */
- if (PseudoZero <= Zero) {
- BadCond(Failure, "Positive expressions can underflow to an\n");
- printf("allegedly negative value\n");
- printf("PseudoZero that prints out as: %g .\n", PseudoZero);
- X = - PseudoZero;
- if (X <= Zero) {
- printf("But -PseudoZero, which should be\n");
- printf("positive, isn't; it prints out as %g .\n", X);
- }
- }
- else {
- BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
- printf("value PseudoZero that prints out as %g .\n", PseudoZero);
- }
- TstPtUf();
- }
- /*=============================================*/
- Milestone = 120;
- /*=============================================*/
- if (CInvrse * Y > CInvrse * Y1) {
- S = H * S;
- E0 = Underflow;
- }
- if (! ((E1 == Zero) || (E1 == E0))) {
- BadCond(Defect, "");
- if (E1 < E0) {
- printf("Products underflow at a higher");
- printf(" threshold than differences.\n");
- if (PseudoZero == Zero)
- E0 = E1;
- }
- else {
- printf("Difference underflows at a higher");
- printf(" threshold than products.\n");
- }
- }
- printf("Smallest strictly positive number found is E0 = %g .\n", E0);
- Z = E0;
- TstPtUf();
- Underflow = E0;
- if (N == 1) Underflow = Y;
- I = 4;
- if (E1 == Zero) I = 3;
- if (UfThold == Zero) I = I - 2;
- UfNGrad = True;
- switch (I) {
- case 1:
- UfThold = Underflow;
- if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
- UfThold = Y;
- BadCond(Failure, "Either accuracy deteriorates as numbers\n");
- printf("approach a threshold = %.17e\n", UfThold);;
- printf(" coming down from %.17e\n", C);
- printf(" or else multiplication gets too many last digits wrong.\n");
- }
- Pause();
- break;
-
- case 2:
- BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
- printf("Q == Y while denying that |Q - Y| == 0; these values\n");
- printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
- printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
- UfThold = Q;
- break;
-
- case 3:
- X = X;
- break;
-
- case 4:
- if ((Q == UfThold) && (E1 == E0)
- && (FABS( UfThold - E1 / E9) <= E1)) {
- UfNGrad = False;
- printf("Underflow is gradual; it incurs Absolute Error =\n");
- printf("(roundoff in UfThold) < E0.\n");
- Y = E0 * CInvrse;
- Y = Y * (OneAndHalf + U2);
- X = CInvrse * (One + U2);
- Y = Y / X;
- IEEE = (Y == E0);
- }
- }
- if (UfNGrad) {
- printf("\n");
- sigsave = sigfpe;
- if (setjmp(ovfl_buf)) {
- printf("Underflow / UfThold failed!\n");
- R = H + H;
- }
- else R = SQRT(Underflow / UfThold);
- sigsave = 0;
- if (R <= H) {
- Z = R * UfThold;
- X = Z * (One + R * H * (One + H));
- }
- else {
- Z = UfThold;
- X = Z * (One + H * H * (One + H));
- }
- if (! ((X == Z) || (X - Z != Zero))) {
- BadCond(Flaw, "");
- printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
- Z9 = X - Z;
- printf("yet X - Z yields %.17e .\n", Z9);
- printf(" Should this NOT signal Underflow, ");
- printf("this is a SERIOUS DEFECT\nthat causes ");
- printf("confusion when innocent statements like\n");;
- printf(" if (X == Z) ... else");
- printf(" ... (f(X) - f(Z)) / (X - Z) ...\n");
- printf("encounter Division by Zero although actually\n");
- sigsave = sigfpe;
- if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
- else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
- sigsave = 0;
- }
- }
- printf("The Underflow threshold is %.17e, %s\n", UfThold,
- " below which");
- printf("calculation may suffer larger Relative error than ");
- printf("merely roundoff.\n");
- Y2 = U1 * U1;
- Y = Y2 * Y2;
- Y2 = Y * U1;
- if (Y2 <= UfThold) {
- if (Y > E0) {
- BadCond(Defect, "");
- I = 5;
- }
- else {
- BadCond(Serious, "");
- I = 4;
- }
- printf("Range is too narrow; U1^%d Underflows.\n", I);
- }
- /*=============================================*/
- /*SPLIT
- }
-#include "paranoia.h"
-part7(){
-*/
- Milestone = 130;
- /*=============================================*/
- Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
- Y2 = Y + Y;
- printf("Since underflow occurs below the threshold\n");
- printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
- printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
- V9 = POW(HInvrse, Y2);
- printf("actually calculating yields: %.17e .\n", V9);
- if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
- BadCond(Serious, "this is not between 0 and underflow\n");
- printf(" threshold = %.17e .\n", UfThold);
- }
- else if (! (V9 > UfThold * (One + E9)))
- printf("This computed value is O.K.\n");
- else {
- BadCond(Defect, "this is not between 0 and underflow\n");
- printf(" threshold = %.17e .\n", UfThold);
- }
- /*=============================================*/
- Milestone = 140;
- /*=============================================*/
- printf("\n");
- /* ...calculate Exp2 == exp(2) == 7.389056099... */
- X = Zero;
- I = 2;
- Y = Two * Three;
- Q = Zero;
- N = 0;
- do {
- Z = X;
- I = I + 1;
- Y = Y / (I + I);
- R = Y + Q;
- X = Z + R;
- Q = (Z - X) + R;
- } while(X > Z);
- Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
- X = Z * Z;
- Exp2 = X * X;
- X = F9;
- Y = X - U1;
- printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
- Exp2);
- for(I = 1;;) {
- Z = X - BInvrse;
- Z = (X + One) / (Z - (One - BInvrse));
- Q = POW(X, Z) - Exp2;
- if (FABS(Q) > TwoForty * U2) {
- N = 1;
- V9 = (X - BInvrse) - (One - BInvrse);
- BadCond(Defect, "Calculated");
- printf(" %.17e for\n", POW(X,Z));
- printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
- printf("\tdiffers from correct value by %.17e .\n", Q);
- printf("\tThis much error may spoil financial\n");
- printf("\tcalculations involving tiny interest rates.\n");
- break;
- }
- else {
- Z = (Y - X) * Two + Y;
- X = Y;
- Y = Z;
- Z = One + (X - F9)*(X - F9);
- if (Z > One && I < NoTrials) I++;
- else {
- if (X > One) {
- if (N == 0)
- printf("Accuracy seems adequate.\n");
- break;
- }
- else {
- X = One + U2;
- Y = U2 + U2;
- Y += X;
- I = 1;
- }
- }
- }
- }
- /*=============================================*/
- Milestone = 150;
- /*=============================================*/
- printf("Testing powers Z^Q at four nearly extreme values.\n");
- N = 0;
- Z = A1;
- Q = FLOOR(Half - LOG(C) / LOG(A1));
- Break = False;
- do {
- X = CInvrse;
- Y = POW(Z, Q);
- IsYeqX();
- Q = - Q;
- X = C;
- Y = POW(Z, Q);
- IsYeqX();
- if (Z < One) Break = True;
- else Z = AInvrse;
- } while ( ! (Break));
- PrintIfNPositive();
- if (N == 0) printf(" ... no discrepancies found.\n");
- printf("\n");
-
- /*=============================================*/
- Milestone = 160;
- /*=============================================*/
- Pause();
- printf("Searching for Overflow threshold:\n");
- printf("This may generate an error.\n");
- Y = - CInvrse;
- V9 = HInvrse * Y;
- sigsave = sigfpe;
- if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
- do {
- V = Y;
- Y = V9;
- V9 = HInvrse * Y;
- } while(V9 < Y);
- I = 1;
-overflow:
- sigsave = 0;
- Z = V9;
- printf("Can `Z = -Y' overflow?\n");
- printf("Trying it on Y = %.17e .\n", Y);
- V9 = - Y;
- V0 = V9;
- if (V - Y == V + V0) printf("Seems O.K.\n");
- else {
- printf("finds a ");
- BadCond(Flaw, "-(-Y) differs from Y.\n");
- }
- if (Z != Y) {
- BadCond(Serious, "");
- printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
- }
- if (I) {
- Y = V * (HInvrse * U2 - HInvrse);
- Z = Y + ((One - HInvrse) * U2) * V;
- if (Z < V0) Y = Z;
- if (Y < V0) V = Y;
- if (V0 - V < V0) V = V0;
- }
- else {
- V = Y * (HInvrse * U2 - HInvrse);
- V = V + ((One - HInvrse) * U2) * Y;
- }
- printf("Overflow threshold is V = %.17e .\n", V);
- if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
- else printf("There is no saturation value because the system traps on overflow.\n");
- V9 = V * One;
- printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
- V9 = V / One;
- printf(" nor for V / 1 = %.17e .\n", V9);
- printf("Any overflow signal separating this * from the one\n");
- printf("above is a DEFECT.\n");
- /*=============================================*/
- Milestone = 170;
- /*=============================================*/
- if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
- BadCond(Failure, "Comparisons involving ");
- printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
- V, V0, UfThold);
- }
- /*=============================================*/
- Milestone = 175;
- /*=============================================*/
- printf("\n");
- for(Indx = 1; Indx <= 3; ++Indx) {
- switch (Indx) {
- case 1: Z = UfThold; break;
- case 2: Z = E0; break;
- case 3: Z = PseudoZero; break;
- }
- if (Z != Zero) {
- V9 = SQRT(Z);
- Y = V9 * V9;
- if (Y / (One - Radix * E9) < Z
- || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
- if (V9 > U1) BadCond(Serious, "");
- else BadCond(Defect, "");
- printf("Comparison alleges that what prints as Z = %.17e\n", Z);
- printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
- }
- }
- }
- /*=============================================*/
- Milestone = 180;
- /*=============================================*/
- for(Indx = 1; Indx <= 2; ++Indx) {
- if (Indx == 1) Z = V;
- else Z = V0;
- V9 = SQRT(Z);
- X = (One - Radix * E9) * V9;
- V9 = V9 * X;
- if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
- Y = V9;
- if (X < W) BadCond(Serious, "");
- else BadCond(Defect, "");
- printf("Comparison alleges that Z = %17e\n", Z);
- printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
- }
- }
- /*=============================================*/
- /*SPLIT
- }
-#include "paranoia.h"
-part8(){
-*/
- Milestone = 190;
- /*=============================================*/
- Pause();
- X = UfThold * V;
- Y = Radix * Radix;
- if (X*Y < One || X > Y) {
- if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
- else BadCond(Flaw, "");
-
- printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
- X, "is too far from 1.\n");
- }
- /*=============================================*/
- Milestone = 200;
- /*=============================================*/
- for (Indx = 1; Indx <= 5; ++Indx) {
- X = F9;
- switch (Indx) {
- case 2: X = One + U2; break;
- case 3: X = V; break;
- case 4: X = UfThold; break;
- case 5: X = Radix;
- }
- Y = X;
- sigsave = sigfpe;
- if (setjmp(ovfl_buf))
- printf(" X / X traps when X = %g\n", X);
- else {
- V9 = (Y / X - Half) - Half;
- if (V9 == Zero) continue;
- if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
- else BadCond(Serious, "");
- printf(" X / X differs from 1 when X = %.17e\n", X);
- printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
- }
- sigsave = 0;
- }
- /*=============================================*/
- Milestone = 210;
- /*=============================================*/
- MyZero = Zero;
- printf("\n");
- printf("What message and/or values does Division by Zero produce?\n") ;
-#ifndef NOPAUSE
- printf("This can interupt your program. You can ");
- printf("skip this part if you wish.\n");
- printf("Do you wish to compute 1 / 0? ");
- fflush(stdout);
- read (KEYBOARD, ch, 8);
- if ((ch[0] == 'Y') || (ch[0] == 'y')) {
-#endif
- sigsave = sigfpe;
- printf(" Trying to compute 1 / 0 produces ...");
- if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero);
- sigsave = 0;
-#ifndef NOPAUSE
- }
- else printf("O.K.\n");
- printf("\nDo you wish to compute 0 / 0? ");
- fflush(stdout);
- read (KEYBOARD, ch, 80);
- if ((ch[0] == 'Y') || (ch[0] == 'y')) {
-#endif
- sigsave = sigfpe;
- printf("\n Trying to compute 0 / 0 produces ...");
- if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero);
- sigsave = 0;
-#ifndef NOPAUSE
- }
- else printf("O.K.\n");
-#endif
- /*=============================================*/
- Milestone = 220;
- /*=============================================*/
- Pause();
- printf("\n");
- {
- static char *msg[] = {
- "FAILUREs encountered =",
- "SERIOUS DEFECTs discovered =",
- "DEFECTs discovered =",
- "FLAWs discovered =" };
- int i;
- for(i = 0; i < 4; i++) if (ErrCnt[i])
- printf("The number of %-29s %d.\n",
- msg[i], ErrCnt[i]);
- }
- printf("\n");
- if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
- + ErrCnt[Flaw]) > 0) {
- if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
- Defect] == 0) && (ErrCnt[Flaw] > 0)) {
- printf("The arithmetic diagnosed seems ");
- printf("Satisfactory though flawed.\n");
- }
- if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
- && ( ErrCnt[Defect] > 0)) {
- printf("The arithmetic diagnosed may be Acceptable\n");
- printf("despite inconvenient Defects.\n");
- }
- if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
- printf("The arithmetic diagnosed has ");
- printf("unacceptable Serious Defects.\n");
- }
- if (ErrCnt[Failure] > 0) {
- printf("Potentially fatal FAILURE may have spoiled this");
- printf(" program's subsequent diagnoses.\n");
- }
- }
- else {
- printf("No failures, defects nor flaws have been discovered.\n");
- if (! ((RMult == Rounded) && (RDiv == Rounded)
- && (RAddSub == Rounded) && (RSqrt == Rounded)))
- printf("The arithmetic diagnosed seems Satisfactory.\n");
- else {
- if (StickyBit >= One &&
- (Radix - Two) * (Radix - Nine - One) == Zero) {
- printf("Rounding appears to conform to ");
- printf("the proposed IEEE standard P");
- if ((Radix == Two) &&
- ((Precision - Four * Three * Two) *
- ( Precision - TwentySeven -
- TwentySeven + One) == Zero))
- printf("754");
- else printf("854");
- if (IEEE) printf(".\n");
- else {
- printf(",\nexcept for possibly Double Rounding");
- printf(" during Gradual Underflow.\n");
- }
- }
- printf("The arithmetic diagnosed appears to be Excellent!\n");
- }
- }
- if (fpecount)
- printf("\nA total of %d floating point exceptions were registered.\n",
- fpecount);
- printf("END OF TEST.\n");
- return 0;
- }
-
-/*SPLIT subs.c
-#include "paranoia.h"
-*/
-
-/* Sign */
-
-FLOAT Sign (X)
-FLOAT X;
-{ return X >= 0. ? 1.0 : -1.0; }
-
-/* Pause */
-
-Pause()
-{
-#ifndef NOPAUSE
- char ch[8];
-
- printf("\nTo continue, press RETURN");
- fflush(stdout);
- read(KEYBOARD, ch, 8);
-#endif
- printf("\nDiagnosis resumes after milestone Number %d", Milestone);
- printf(" Page: %d\n\n", PageNo);
- ++Milestone;
- ++PageNo;
- }
-
- /* TstCond */
-
-TstCond (K, Valid, T)
-int K, Valid;
-char *T;
-{ if (! Valid) { BadCond(K,T); printf(".\n"); } }
-
-BadCond(K, T)
-int K;
-char *T;
-{
- static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
-
- ErrCnt [K] = ErrCnt [K] + 1;
- printf("%s: %s", msg[K], T);
- }
-
-/* Random */
-/* Random computes
- X = (Random1 + Random9)^5
- Random1 = X - FLOOR(X) + 0.000005 * X;
- and returns the new value of Random1
-*/
-
-FLOAT Random()
-{
- FLOAT X, Y;
-
- X = Random1 + Random9;
- Y = X * X;
- Y = Y * Y;
- X = X * Y;
- Y = X - FLOOR(X);
- Random1 = Y + X * 0.000005;
- return(Random1);
- }
-
-/* SqXMinX */
-
-SqXMinX (ErrKind)
-int ErrKind;
-{
- FLOAT XA, XB;
-
- XB = X * BInvrse;
- XA = X - XB;
- SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
- if (SqEr != Zero) {
- if (SqEr < MinSqEr) MinSqEr = SqEr;
- if (SqEr > MaxSqEr) MaxSqEr = SqEr;
- J = J + 1.0;
- BadCond(ErrKind, "\n");
- printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr);
- printf("\tinstead of correct value 0 .\n");
- }
- }
-
-/* NewD */
-
-NewD()
-{
- X = Z1 * Q;
- X = FLOOR(Half - X / Radix) * Radix + X;
- Q = (Q - X * Z) / Radix + X * X * (D / Radix);
- Z = Z - Two * X * D;
- if (Z <= Zero) {
- Z = - Z;
- Z1 = - Z1;
- }
- D = Radix * D;
- }
-
-/* SR3750 */
-
-SR3750()
-{
- if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
- I = I + 1;
- X2 = SQRT(X * D);
- Y2 = (X2 - Z2) - (Y - Z2);
- X2 = X8 / (Y - Half);
- X2 = X2 - Half * X2 * X2;
- SqEr = (Y2 + Half) + (Half - X2);
- if (SqEr < MinSqEr) MinSqEr = SqEr;
- SqEr = Y2 - X2;
- if (SqEr > MaxSqEr) MaxSqEr = SqEr;
- }
- }
-
-/* IsYeqX */
-
-IsYeqX()
-{
- if (Y != X) {
- if (N <= 0) {
- if (Z == Zero && Q <= Zero)
- printf("WARNING: computing\n");
- else BadCond(Defect, "computing\n");
- printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
- printf("\tyielded %.17e;\n", Y);
- printf("\twhich compared unequal to correct %.17e ;\n",
- X);
- printf("\t\tthey differ by %.17e .\n", Y - X);
- }
- N = N + 1; /* ... count discrepancies. */
- }
- }
-
-/* SR3980 */
-
-SR3980()
-{
- do {
- Q = (FLOAT) I;
- Y = POW(Z, Q);
- IsYeqX();
- if (++I > M) break;
- X = Z * X;
- } while ( X < W );
- }
-
-/* PrintIfNPositive */
-
-PrintIfNPositive()
-{
- if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
- }
-
-/* TstPtUf */
-
-TstPtUf()
-{
- N = 0;
- if (Z != Zero) {
- printf("Since comparison denies Z = 0, evaluating ");
- printf("(Z + Z) / Z should be safe.\n");
- sigsave = sigfpe;
- if (setjmp(ovfl_buf)) goto very_serious;
- Q9 = (Z + Z) / Z;
- printf("What the machine gets for (Z + Z) / Z is %.17e .\n",
- Q9);
- if (FABS(Q9 - Two) < Radix * U2) {
- printf("This is O.K., provided Over/Underflow");
- printf(" has NOT just been signaled.\n");
- }
- else {
- if ((Q9 < One) || (Q9 > Two)) {
-very_serious:
- N = 1;
- ErrCnt [Serious] = ErrCnt [Serious] + 1;
- printf("This is a VERY SERIOUS DEFECT!\n");
- }
- else {
- N = 1;
- ErrCnt [Defect] = ErrCnt [Defect] + 1;
- printf("This is a DEFECT!\n");
- }
- }
- sigsave = 0;
- V9 = Z * One;
- Random1 = V9;
- V9 = One * Z;
- Random2 = V9;
- V9 = Z / One;
- if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
- if (N > 0) Pause();
- }
- else {
- N = 1;
- BadCond(Defect, "What prints as Z = ");
- printf("%.17e\n\tcompares different from ", Z);
- if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
- if (! ((Z == Random2)
- || (Random2 == Random1)))
- printf("1 * Z == %g\n", Random2);
- if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
- if (Random2 != Random1) {
- ErrCnt [Defect] = ErrCnt [Defect] + 1;
- BadCond(Defect, "Multiplication does not commute!\n");
- printf("\tComparison alleges that 1 * Z = %.17e\n",
- Random2);
- printf("\tdiffers from Z * 1 = %.17e\n", Random1);
- }
- Pause();
- }
- }
- }
-
-notify(s)
-char *s;
-{
- printf("%s test appears to be inconsistent...\n", s);
- printf(" PLEASE NOTIFY KARPINKSI!\n");
- }
-
-/*SPLIT msgs.c */
-
-/* Instructions */
-
-msglist(s)
-char **s;
-{ while(*s) printf("%s\n", *s++); }
-
-Instructions()
-{
- static char *instr[] = {
- "Lest this program stop prematurely, i.e. before displaying\n",
- " `END OF TEST',\n",
- "try to persuade the computer NOT to terminate execution when an",
- "error like Over/Underflow or Division by Zero occurs, but rather",
- "to persevere with a surrogate value after, perhaps, displaying some",
- "warning. If persuasion avails naught, don't despair but run this",
- "program anyway to see how many milestones it passes, and then",
- "amend it to make further progress.\n",
- "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
- 0};
-
- msglist(instr);
- }
-
-/* Heading */
-
-Heading()
-{
- static char *head[] = {
- "Users are invited to help debug and augment this program so it will",
- "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
- "Please send suggestions and interesting results to",
- "\tRichard Karpinski",
- "\tComputer Center U-76",
- "\tUniversity of California",
- "\tSan Francisco, CA 94143-0704, USA\n",
- "In doing so, please include the following information:",
-#ifdef Single
- "\tPrecision:\tsingle;",
-#else
- "\tPrecision:\tdouble;",
-#endif
- "\tVersion:\t10 February 1989;",
- "\tComputer:\n",
- "\tCompiler:\n",
- "\tOptimization level:\n",
- "\tOther relevant compiler options:",
- 0};
-
- msglist(head);
- }
-
-/* Characteristics */
-
-Characteristics()
-{
- static char *chars[] = {
- "Running this program should reveal these characteristics:",
- " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
- " Precision = number of significant digits carried.",
- " U2 = Radix/Radix^Precision = One Ulp",
- "\t(OneUlpnit in the Last Place) of 1.000xxx .",
- " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
- " Adequacy of guard digits for Mult., Div. and Subt.",
- " Whether arithmetic is chopped, correctly rounded, or something else",
- "\tfor Mult., Div., Add/Subt. and Sqrt.",
- " Whether a Sticky Bit used correctly for rounding.",
- " UnderflowThreshold = an underflow threshold.",
- " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
- " V = an overflow threshold, roughly.",
- " V0 tells, roughly, whether Infinity is represented.",
- " Comparisions are checked for consistency with subtraction",
- "\tand for contamination with pseudo-zeros.",
- " Sqrt is tested. Y^X is not tested.",
- " Extra-precise subexpressions are revealed but NOT YET tested.",
- " Decimal-Binary conversion is NOT YET tested for accuracy.",
- 0};
-
- msglist(chars);
- }
-
-History()
-
-{ /* History */
- /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
- with further massaging by David M. Gay. */
-
- static char *hist[] = {
- "The program attempts to discriminate among",
- " FLAWs, like lack of a sticky bit,",
- " Serious DEFECTs, like lack of a guard digit, and",
- " FAILUREs, like 2+2 == 5 .",
- "Failures may confound subsequent diagnoses.\n",
- "The diagnostic capabilities of this program go beyond an earlier",
- "program called `MACHAR', which can be found at the end of the",
- "book `Software Manual for the Elementary Functions' (1980) by",
- "W. J. Cody and W. Waite. Although both programs try to discover",
- "the Radix, Precision and range (over/underflow thresholds)",
- "of the arithmetic, this program tries to cope with a wider variety",
- "of pathologies, and to say how well the arithmetic is implemented.",
- "\nThe program is based upon a conventional radix representation for",
- "floating-point numbers, but also allows logarithmic encoding",
- "as used by certain early WANG machines.\n",
- "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
- "see source comments for more history.",
- 0};
-
- msglist(hist);
- }
-
-double
-pow(x, y) /* return x ^ y (exponentiation) */
-double x, y;
-{
- extern double exp(), frexp(), ldexp(), log(), modf();
- double xy, ye;
- long i;
- int ex, ey = 0, flip = 0;
-
- if (!y) return 1.0;
-
- if ((y < -1100. || y > 1100.) && x != -1.) return exp(y * log(x));
-
- if (y < 0.) { y = -y; flip = 1; }
- y = modf(y, &ye);
- if (y) xy = exp(y * log(x));
- else xy = 1.0;
- /* next several lines assume >= 32 bit integers */
- x = frexp(x, &ex);
- if (i = ye) for(;;) {
- if (i & 1) { xy *= x; ey += ex; }
- if (!(i >>= 1)) break;
- x *= x;
- ex *= 2;
- if (x < .5) { x *= 2.; ex -= 1; }
- }
- if (flip) { xy = 1. / xy; ey = -ey; }
- return ldexp(xy, ey);
-}
+#undef V9 +#define NOPAUSE +/* A C version of Kahan's Floating Point Test "Paranoia" + + Thos Sumner, UCSF, Feb. 1985 + David Gay, BTL, Jan. 1986 + + This is a rewrite from the Pascal version by + + B. A. Wichmann, 18 Jan. 1985 + + (and does NOT exhibit good C programming style). + +(C) Apr 19 1983 in BASIC version by: + Professor W. M. Kahan, + 567 Evans Hall + Electrical Engineering & Computer Science Dept. + University of California + Berkeley, California 94720 + USA + +converted to Pascal by: + B. A. Wichmann + National Physical Laboratory + Teddington Middx + TW11 OLW + UK + +converted to C by: + + David M. Gay and Thos Sumner + AT&T Bell Labs Computer Center, Rm. U-76 + 600 Mountain Avenue University of California + Murray Hill, NJ 07974 San Francisco, CA 94143 + USA USA + +with simultaneous corrections to the Pascal source (reflected +in the Pascal source available over netlib). +[A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.] + +Reports of results on various systems from all the versions +of Paranoia are being collected by Richard Karpinski at the +same address as Thos Sumner. This includes sample outputs, +bug reports, and criticisms. + +You may copy this program freely if you acknowledge its source. +Comments on the Pascal version to NPL, please. + + +The C version catches signals from floating-point exceptions. +If signal(SIGFPE,...) is unavailable in your environment, you may +#define NOSIGNAL to comment out the invocations of signal. + +This source file is too big for some C compilers, but may be split +into pieces. Comments containing "SPLIT" suggest convenient places +for this splitting. At the end of these comments is an "ed script" +(for the UNIX(tm) editor ed) that will do this splitting. + +By #defining Single when you compile this source, you may obtain +a single-precision C version of Paranoia. + + +The following is from the introductory commentary from Wichmann's work: + +The BASIC program of Kahan is written in Microsoft BASIC using many +facilities which have no exact analogy in Pascal. The Pascal +version below cannot therefore be exactly the same. Rather than be +a minimal transcription of the BASIC program, the Pascal coding +follows the conventional style of block-structured languages. Hence +the Pascal version could be useful in producing versions in other +structured languages. + +Rather than use identifiers of minimal length (which therefore have +little mnemonic significance), the Pascal version uses meaningful +identifiers as follows [Note: A few changes have been made for C]: + + +BASIC C BASIC C BASIC C + + A J S StickyBit + A1 AInverse J0 NoErrors T + B Radix [Failure] T0 Underflow + B1 BInverse J1 NoErrors T2 ThirtyTwo + B2 RadixD2 [SeriousDefect] T5 OneAndHalf + B9 BMinusU2 J2 NoErrors T7 TwentySeven + C [Defect] T8 TwoForty + C1 CInverse J3 NoErrors U OneUlp + D [Flaw] U0 UnderflowThreshold + D4 FourD K PageNo U1 + E0 L Milestone U2 + E1 M V + E2 Exp2 N V0 + E3 N1 V8 + E5 MinSqEr O Zero V9 + E6 SqEr O1 One W + E7 MaxSqEr O2 Two X + E8 O3 Three X1 + E9 O4 Four X8 + F1 MinusOne O5 Five X9 Random1 + F2 Half O8 Eight Y + F3 Third O9 Nine Y1 + F6 P Precision Y2 + F9 Q Y9 Random2 + G1 GMult Q8 Z + G2 GDiv Q9 Z0 PseudoZero + G3 GAddSub R Z1 + H R1 RMult Z2 + H1 HInverse R2 RDiv Z9 + I R3 RAddSub + IO NoTrials R4 RSqrt + I3 IEEE R9 Random9 + + SqRWrng + +All the variables in BASIC are true variables and in consequence, +the program is more difficult to follow since the "constants" must +be determined (the glossary is very helpful). The Pascal version +uses Real constants, but checks are added to ensure that the values +are correctly converted by the compiler. + +The major textual change to the Pascal version apart from the +identifiersis that named procedures are used, inserting parameters +wherehelpful. New procedures are also introduced. The +correspondence is as follows: + + +BASIC Pascal +lines + + 90- 140 Pause + 170- 250 Instructions + 380- 460 Heading + 480- 670 Characteristics + 690- 870 History +2940-2950 Random +3710-3740 NewD +4040-4080 DoesYequalX +4090-4110 PrintIfNPositive +4640-4850 TestPartialUnderflow + +=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*= + +Below is an "ed script" that splits para.c into 10 files +of the form part[1-8].c, subs.c, and msgs.c, plus a header +file, paranoia.h, that these files require. + +r paranoia.c +$ +?SPLIT ++,$w msgs.c + .,$d +?SPLIT + .d ++d +-,$w subs.c +-,$d +?part8 ++d +?include + .,$w part8.c + .,$d +-d +?part7 ++d +?include + .,$w part7.c + .,$d +-d +?part6 ++d +?include + .,$w part6.c + .,$d +-d +?part5 ++d +?include + .,$w part5.c + .,$d +-d +?part4 ++d +?include + .,$w part4.c + .,$d +-d +?part3 ++d +?include + .,$w part3.c + .,$d +-d +?part2 ++d +?include + .,$w part2.c + .,$d +?SPLIT + .d +1,/^#include/-1d +1,$w part1.c +/Computed constants/,$d +1,$s/^int/extern &/ +1,$s/^FLOAT/extern &/ +1,$s/^char/extern &/ +1,$s! = .*!;! +/^Guard/,/^Round/s/^/extern / +/^jmp_buf/s/^/extern / +/^Sig_type/s/^/extern / +s/$/\ +extern void sigfpe();/ +w paranoia.h +q + +*/ + +#include <stdio.h> +#ifndef NOSIGNAL +#include <signal.h> +#endif +#include <setjmp.h> + +extern double fabs(), floor(), log(), pow(), sqrt(); + +#ifdef Single +#define FLOAT float +#define FABS(x) (float)fabs((double)(x)) +#define FLOOR(x) (float)floor((double)(x)) +#define LOG(x) (float)log((double)(x)) +#define POW(x,y) (float)pow((double)(x),(double)(y)) +#define SQRT(x) (float)sqrt((double)(x)) +#else +#define FLOAT double +#define FABS(x) fabs(x) +#define FLOOR(x) floor(x) +#define LOG(x) log(x) +#define POW(x,y) pow(x,y) +#define SQRT(x) sqrt(x) +#endif + +jmp_buf ovfl_buf; +typedef void (*Sig_type)(); +Sig_type sigsave; + +#define KEYBOARD 0 + +FLOAT Radix, BInvrse, RadixD2, BMinusU2; +FLOAT Sign(), Random(); + +/*Small floating point constants.*/ +FLOAT Zero = 0.0; +FLOAT Half = 0.5; +FLOAT One = 1.0; +FLOAT Two = 2.0; +FLOAT Three = 3.0; +FLOAT Four = 4.0; +FLOAT Five = 5.0; +FLOAT Eight = 8.0; +FLOAT Nine = 9.0; +FLOAT TwentySeven = 27.0; +FLOAT ThirtyTwo = 32.0; +FLOAT TwoForty = 240.0; +FLOAT MinusOne = -1.0; +FLOAT OneAndHalf = 1.5; +/*Integer constants*/ +int NoTrials = 20; /*Number of tests for commutativity. */ +#define False 0 +#define True 1 + +/* Definitions for declared types + Guard == (Yes, No); + Rounding == (Chopped, Rounded, Other); + Message == packed array [1..40] of char; + Class == (Flaw, Defect, Serious, Failure); + */ +#define Yes 1 +#define No 0 +#define Chopped 2 +#define Rounded 1 +#define Other 0 +#define Flaw 3 +#define Defect 2 +#define Serious 1 +#define Failure 0 +typedef int Guard, Rounding, Class; +typedef char Message; + +/* Declarations of Variables */ +int Indx; +char ch[8]; +FLOAT AInvrse, A1; +FLOAT C, CInvrse; +FLOAT D, FourD; +FLOAT E0, E1, Exp2, E3, MinSqEr; +FLOAT SqEr, MaxSqEr, E9; +FLOAT Third; +FLOAT F6, F9; +FLOAT H, HInvrse; +int I; +FLOAT StickyBit, J; +FLOAT MyZero; +FLOAT Precision; +FLOAT Q, Q9; +FLOAT R, Random9; +FLOAT T, Underflow, S; +FLOAT OneUlp, UfThold, U1, U2; +FLOAT V, V0, V9; +FLOAT W; +FLOAT X, X1, X2, X8, Random1; +FLOAT Y, Y1, Y2, Random2; +FLOAT Z, PseudoZero, Z1, Z2, Z9; +int ErrCnt[4]; +int fpecount; +int Milestone; +int PageNo; +int M, N, N1; +Guard GMult, GDiv, GAddSub; +Rounding RMult, RDiv, RAddSub, RSqrt; +int Break, Done, NotMonot, Monot, Anomaly, IEEE, + SqRWrng, UfNGrad; +/* Computed constants. */ +/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */ +/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */ + +/* floating point exception receiver */ + void +sigfpe(i) +{ + fpecount++; + printf("\n* * * FLOATING-POINT ERROR * * *\n"); + fflush(stdout); + if (sigsave) { +#ifndef NOSIGNAL + signal(SIGFPE, sigsave); +#endif + sigsave = 0; + longjmp(ovfl_buf, 1); + } + abort(); +} + +main() +{ +#ifdef mc + char *out; + ieee_flags("set", "precision", "double", &out); +#endif + /* First two assignments use integer right-hand sides. */ + Zero = 0; + One = 1; + Two = One + One; + Three = Two + One; + Four = Three + One; + Five = Four + One; + Eight = Four + Four; + Nine = Three * Three; + TwentySeven = Nine * Three; + ThirtyTwo = Four * Eight; + TwoForty = Four * Five * Three * Four; + MinusOne = -One; + Half = One / Two; + OneAndHalf = One + Half; + ErrCnt[Failure] = 0; + ErrCnt[Serious] = 0; + ErrCnt[Defect] = 0; + ErrCnt[Flaw] = 0; + PageNo = 1; + /*=============================================*/ + Milestone = 0; + /*=============================================*/ +#ifndef NOSIGNAL + signal(SIGFPE, sigfpe); +#endif + Instructions(); + Pause(); + Heading(); + Pause(); + Characteristics(); + Pause(); + History(); + Pause(); + /*=============================================*/ + Milestone = 7; + /*=============================================*/ + printf("Program is now RUNNING tests on small integers:\n"); + + TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero) + && (One > Zero) && (One + One == Two), + "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2"); + Z = - Zero; + if (Z != 0.0) { + ErrCnt[Failure] = ErrCnt[Failure] + 1; + printf("Comparison alleges that -0.0 is Non-zero!\n"); + U1 = 0.001; + Radix = 1; + TstPtUf(); + } + TstCond (Failure, (Three == Two + One) && (Four == Three + One) + && (Four + Two * (- Two) == Zero) + && (Four - Three - One == Zero), + "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0"); + TstCond (Failure, (MinusOne == (0 - One)) + && (MinusOne + One == Zero ) && (One + MinusOne == Zero) + && (MinusOne + FABS(One) == Zero) + && (MinusOne + MinusOne * MinusOne == Zero), + "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0"); + TstCond (Failure, Half + MinusOne + Half == Zero, + "1/2 + (-1) + 1/2 != 0"); + /*=============================================*/ + /*SPLIT + part2(); + part3(); + part4(); + part5(); + part6(); + part7(); + part8(); + } +#include "paranoia.h" +part2(){ +*/ + Milestone = 10; + /*=============================================*/ + TstCond (Failure, (Nine == Three * Three) + && (TwentySeven == Nine * Three) && (Eight == Four + Four) + && (ThirtyTwo == Eight * Four) + && (ThirtyTwo - TwentySeven - Four - One == Zero), + "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0"); + TstCond (Failure, (Five == Four + One) && + (TwoForty == Four * Five * Three * Four) + && (TwoForty / Three - Four * Four * Five == Zero) + && ( TwoForty / Four - Five * Three * Four == Zero) + && ( TwoForty / Five - Four * Three * Four == Zero), + "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48"); + if (ErrCnt[Failure] == 0) { + printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n"); + printf("\n"); + } + printf("Searching for Radix and Precision.\n"); + W = One; + do { + W = W + W; + Y = W + One; + Z = Y - W; + Y = Z - One; + } while (MinusOne + FABS(Y) < Zero); + /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/ + Precision = Zero; + Y = One; + do { + Radix = W + Y; + Y = Y + Y; + Radix = Radix - W; + } while ( Radix == Zero); + if (Radix < Two) Radix = One; + printf("Radix = %f .\n", Radix); + if (Radix != 1) { + W = One; + do { + Precision = Precision + One; + W = W * Radix; + Y = W + One; + } while ((Y - W) == One); + } + /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 + ...*/ + U1 = One / W; + U2 = Radix * U1; + printf("Closest relative separation found is U1 = %.7e .\n\n", U1); + printf("Recalculating radix and precision\n "); + + /*save old values*/ + E0 = Radix; + E1 = U1; + E9 = U2; + E3 = Precision; + + X = Four / Three; + Third = X - One; + F6 = Half - Third; + X = F6 + F6; + X = FABS(X - Third); + if (X < U2) X = U2; + + /*... now X = (unknown no.) ulps of 1+...*/ + do { + U2 = X; + Y = Half * U2 + ThirtyTwo * U2 * U2; + Y = One + Y; + X = Y - One; + } while ( ! ((U2 <= X) || (X <= Zero))); + + /*... now U2 == 1 ulp of 1 + ... */ + X = Two / Three; + F6 = X - Half; + Third = F6 + F6; + X = Third - Half; + X = FABS(X + F6); + if (X < U1) X = U1; + + /*... now X == (unknown no.) ulps of 1 -... */ + do { + U1 = X; + Y = Half * U1 + ThirtyTwo * U1 * U1; + Y = Half - Y; + X = Half + Y; + Y = Half - X; + X = Half + Y; + } while ( ! ((U1 <= X) || (X <= Zero))); + /*... now U1 == 1 ulp of 1 - ... */ + if (U1 == E1) printf("confirms closest relative separation U1 .\n"); + else printf("gets better closest relative separation U1 = %.7e .\n", U1); + W = One / U1; + F9 = (Half - U1) + Half; + Radix = FLOOR(0.01 + U2 / U1); + if (Radix == E0) printf("Radix confirmed.\n"); + else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix); + TstCond (Defect, Radix <= Eight + Eight, + "Radix is too big: roundoff problems"); + TstCond (Flaw, (Radix == Two) || (Radix == 10) + || (Radix == One), "Radix is not as good as 2 or 10"); + /*=============================================*/ + Milestone = 20; + /*=============================================*/ + TstCond (Failure, F9 - Half < Half, + "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?"); + X = F9; + I = 1; + Y = X - Half; + Z = Y - Half; + TstCond (Failure, (X != One) + || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0"); + X = One + U2; + I = 0; + /*=============================================*/ + Milestone = 25; + /*=============================================*/ + /*... BMinusU2 = nextafter(Radix, 0) */ + BMinusU2 = Radix - One; + BMinusU2 = (BMinusU2 - U2) + One; + /* Purify Integers */ + if (Radix != One) { + X = - TwoForty * LOG(U1) / LOG(Radix); + Y = FLOOR(Half + X); + if (FABS(X - Y) * Four < One) X = Y; + Precision = X / TwoForty; + Y = FLOOR(Half + Precision); + if (FABS(Precision - Y) * TwoForty < Half) Precision = Y; + } + if ((Precision != FLOOR(Precision)) || (Radix == One)) { + printf("Precision cannot be characterized by an Integer number\n"); + printf("of significant digits but, by itself, this is a minor flaw.\n"); + } + if (Radix == One) + printf("logarithmic encoding has precision characterized solely by U1.\n"); + else printf("The number of significant digits of the Radix is %f .\n", + Precision); + TstCond (Serious, U2 * Nine * Nine * TwoForty < One, + "Precision worse than 5 decimal figures "); + /*=============================================*/ + Milestone = 30; + /*=============================================*/ + /* Test for extra-precise subepressions */ + X = FABS(((Four / Three - One) - One / Four) * Three - One / Four); + do { + Z2 = X; + X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One; + } while ( ! ((Z2 <= X) || (X <= Zero))); + X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four); + do { + Z1 = Z; + Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1)) + + One / Two)) + One / Two; + } while ( ! ((Z1 <= Z) || (Z <= Zero))); + do { + do { + Y1 = Y; + Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half + )) + Half; + } while ( ! ((Y1 <= Y) || (Y <= Zero))); + X1 = X; + X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9; + } while ( ! ((X1 <= X) || (X <= Zero))); + if ((X1 != Y1) || (X1 != Z1)) { + BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n"); + printf("respectively %.7e, %.7e, %.7e,\n", X1, Y1, Z1); + printf("are symptoms of inconsistencies introduced\n"); + printf("by extra-precise evaluation of arithmetic subexpressions.\n"); + notify("Possibly some part of this"); + if ((X1 == U1) || (Y1 == U1) || (Z1 == U1)) printf( + "That feature is not tested further by this program.\n") ; + } + else { + if ((Z1 != U1) || (Z2 != U2)) { + if ((Z1 >= U1) || (Z2 >= U2)) { + BadCond(Failure, ""); + notify("Precision"); + printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1); + printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2); + } + else { + if ((Z1 <= Zero) || (Z2 <= Zero)) { + printf("Because of unusual Radix = %f", Radix); + printf(", or exact rational arithmetic a result\n"); + printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2); + notify("of an\nextra-precision"); + } + if (Z1 != Z2 || Z1 > Zero) { + X = Z1 / U1; + Y = Z2 / U2; + if (Y > X) X = Y; + Q = - LOG(X); + printf("Some subexpressions appear to be calculated extra\n"); + printf("precisely with about %g extra B-digits, i.e.\n", + (Q / LOG(Radix))); + printf("roughly %g extra significant decimals.\n", + Q / LOG(10.)); + } + printf("That feature is not tested further by this program.\n"); + } + } + } + Pause(); + /*=============================================*/ + /*SPLIT + } +#include "paranoia.h" +part3(){ +*/ + Milestone = 35; + /*=============================================*/ + if (Radix >= Two) { + X = W / (Radix * Radix); + Y = X + One; + Z = Y - X; + T = Z + U2; + X = T - Z; + TstCond (Failure, X == U2, + "Subtraction is not normalized X=Y,X+Z != Y+Z!"); + if (X == U2) printf( + "Subtraction appears to be normalized, as it should be."); + } + printf("\nChecking for guard digit in *, /, and -.\n"); + Y = F9 * One; + Z = One * F9; + X = F9 - Half; + Y = (Y - Half) - X; + Z = (Z - Half) - X; + X = One + U2; + T = X * Radix; + R = Radix * X; + X = T - Radix; + X = X - Radix * U2; + T = R - Radix; + T = T - Radix * U2; + X = X * (Radix - One); + T = T * (Radix - One); + if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes; + else { + GMult = No; + TstCond (Serious, False, + "* lacks a Guard Digit, so 1*X != X"); + } + Z = Radix * U2; + X = One + Z; + Y = FABS((X + Z) - X * X) - U2; + X = One - U2; + Z = FABS((X - U2) - X * X) - U1; + TstCond (Failure, (Y <= Zero) + && (Z <= Zero), "* gets too many final digits wrong.\n"); + Y = One - U2; + X = One + U2; + Z = One / Y; + Y = Z - X; + X = One / Three; + Z = Three / Nine; + X = X - Z; + T = Nine / TwentySeven; + Z = Z - T; + TstCond(Defect, X == Zero && Y == Zero && Z == Zero, + "Division lacks a Guard Digit, so error can exceed 1 ulp\nor 1/3 and 3/9 and 9/27 may disagree"); + Y = F9 / One; + X = F9 - Half; + Y = (Y - Half) - X; + X = One + U2; + T = X / One; + X = T - X; + if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes; + else { + GDiv = No; + TstCond (Serious, False, + "Division lacks a Guard Digit, so X/1 != X"); + } + X = One / (One + U2); + Y = X - Half - Half; + TstCond (Serious, Y < Zero, + "Computed value of 1/1.000..1 >= 1"); + X = One - U2; + Y = One + Radix * U2; + Z = X * Radix; + T = Y * Radix; + R = Z / Radix; + StickyBit = T / Radix; + X = R - X; + Y = StickyBit - Y; + TstCond (Failure, X == Zero && Y == Zero, + "* and/or / gets too many last digits wrong"); + Y = One - U1; + X = One - F9; + Y = One - Y; + T = Radix - U2; + Z = Radix - BMinusU2; + T = Radix - T; + if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes; + else { + GAddSub = No; + TstCond (Serious, False, + "- lacks Guard Digit, so cancellation is obscured"); + } + if (F9 != One && F9 - One >= Zero) { + BadCond(Serious, "comparison alleges (1-U1) < 1 although\n"); + printf(" subtraction yields (1-U1) - 1 = 0 , thereby vitiating\n"); + printf(" such precautions against division by zero as\n"); + printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n"); + } + if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf( + " *, /, and - appear to have guard digits, as they should.\n"); + /*=============================================*/ + Milestone = 40; + /*=============================================*/ + Pause(); + printf("Checking rounding on multiply, divide and add/subtract.\n"); + RMult = Other; + RDiv = Other; + RAddSub = Other; + RadixD2 = Radix / Two; + A1 = Two; + Done = False; + do { + AInvrse = Radix; + do { + X = AInvrse; + AInvrse = AInvrse / A1; + } while ( ! (FLOOR(AInvrse) != AInvrse)); + Done = (X == One) || (A1 > Three); + if (! Done) A1 = Nine + One; + } while ( ! (Done)); + if (X == One) A1 = Radix; + AInvrse = One / A1; + X = A1; + Y = AInvrse; + Done = False; + do { + Z = X * Y - Half; + TstCond (Failure, Z == Half, + "X * (1/X) differs from 1"); + Done = X == Radix; + X = Radix; + Y = One / X; + } while ( ! (Done)); + Y2 = One + U2; + Y1 = One - U2; + X = OneAndHalf - U2; + Y = OneAndHalf + U2; + Z = (X - U2) * Y2; + T = Y * Y1; + Z = Z - X; + T = T - X; + X = X * Y2; + Y = (Y + U2) * Y1; + X = X - OneAndHalf; + Y = Y - OneAndHalf; + if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) { + X = (OneAndHalf + U2) * Y2; + Y = OneAndHalf - U2 - U2; + Z = OneAndHalf + U2 + U2; + T = (OneAndHalf - U2) * Y1; + X = X - (Z + U2); + StickyBit = Y * Y1; + S = Z * Y2; + T = T - Y; + Y = (U2 - Y) + StickyBit; + Z = S - (Z + U2 + U2); + StickyBit = (Y2 + U2) * Y1; + Y1 = Y2 * Y1; + StickyBit = StickyBit - Y2; + Y1 = Y1 - Half; + if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero) + && ( StickyBit == Zero) && (Y1 == Half)) { + RMult = Rounded; + printf("Multiplication appears to round correctly.\n"); + } + else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero) + && (T < Zero) && (StickyBit + U2 == Zero) + && (Y1 < Half)) { + RMult = Chopped; + printf("Multiplication appears to chop.\n"); + } + else printf("* is neither chopped nor correctly rounded.\n"); + if ((RMult == Rounded) && (GMult == No)) notify("Multiplication"); + } + else printf("* is neither chopped nor correctly rounded.\n"); + /*=============================================*/ + Milestone = 45; + /*=============================================*/ + Y2 = One + U2; + Y1 = One - U2; + Z = OneAndHalf + U2 + U2; + X = Z / Y2; + T = OneAndHalf - U2 - U2; + Y = (T - U2) / Y1; + Z = (Z + U2) / Y2; + X = X - OneAndHalf; + Y = Y - T; + T = T / Y1; + Z = Z - (OneAndHalf + U2); + T = (U2 - OneAndHalf) + T; + if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) { + X = OneAndHalf / Y2; + Y = OneAndHalf - U2; + Z = OneAndHalf + U2; + X = X - Y; + T = OneAndHalf / Y1; + Y = Y / Y1; + T = T - (Z + U2); + Y = Y - Z; + Z = Z / Y2; + Y1 = (Y2 + U2) / Y2; + Z = Z - OneAndHalf; + Y2 = Y1 - Y2; + Y1 = (F9 - U1) / F9; + if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero) + && (Y2 == Zero) && (Y2 == Zero) + && (Y1 - Half == F9 - Half )) { + RDiv = Rounded; + printf("Division appears to round correctly.\n"); + if (GDiv == No) notify("Division"); + } + else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero) + && (Y2 < Zero) && (Y1 - Half < F9 - Half)) { + RDiv = Chopped; + printf("Division appears to chop.\n"); + } + } + if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n"); + BInvrse = One / Radix; + TstCond (Failure, (BInvrse * Radix - Half == Half), + "Radix * ( 1 / Radix ) differs from 1"); + /*=============================================*/ + /*SPLIT + } +#include "paranoia.h" +part4(){ +*/ + Milestone = 50; + /*=============================================*/ + TstCond (Failure, ((F9 + U1) - Half == Half) + && ((BMinusU2 + U2 ) - One == Radix - One), + "Incomplete carry-propagation in Addition"); + X = One - U1 * U1; + Y = One + U2 * (One - U2); + Z = F9 - Half; + X = (X - Half) - Z; + Y = Y - One; + if ((X == Zero) && (Y == Zero)) { + RAddSub = Chopped; + printf("Add/Subtract appears to be chopped.\n"); + } + if (GAddSub == Yes) { + X = (Half + U2) * U2; + Y = (Half - U2) * U2; + X = One + X; + Y = One + Y; + X = (One + U2) - X; + Y = One - Y; + if ((X == Zero) && (Y == Zero)) { + X = (Half + U2) * U1; + Y = (Half - U2) * U1; + X = One - X; + Y = One - Y; + X = F9 - X; + Y = One - Y; + if ((X == Zero) && (Y == Zero)) { + RAddSub = Rounded; + printf("Addition/Subtraction appears to round correctly.\n"); + if (GAddSub == No) notify("Add/Subtract"); + } + else printf("Addition/Subtraction neither rounds nor chops.\n"); + } + else printf("Addition/Subtraction neither rounds nor chops.\n"); + } + else printf("Addition/Subtraction neither rounds nor chops.\n"); + S = One; + X = One + Half * (One + Half); + Y = (One + U2) * Half; + Z = X - Y; + T = Y - X; + StickyBit = Z + T; + if (StickyBit != Zero) { + S = Zero; + BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n"); + } + StickyBit = Zero; + if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes) + && (RMult == Rounded) && (RDiv == Rounded) + && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) { + printf("Checking for sticky bit.\n"); + X = (Half + U1) * U2; + Y = Half * U2; + Z = One + Y; + T = One + X; + if ((Z - One <= Zero) && (T - One >= U2)) { + Z = T + Y; + Y = Z - X; + if ((Z - T >= U2) && (Y - T == Zero)) { + X = (Half + U1) * U1; + Y = Half * U1; + Z = One - Y; + T = One - X; + if ((Z - One == Zero) && (T - F9 == Zero)) { + Z = (Half - U1) * U1; + T = F9 - Z; + Q = F9 - Y; + if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) { + Z = (One + U2) * OneAndHalf; + T = (OneAndHalf + U2) - Z + U2; + X = One + Half / Radix; + Y = One + Radix * U2; + Z = X * Y; + if (T == Zero && X + Radix * U2 - Z == Zero) { + if (Radix != Two) { + X = Two + U2; + Y = X / Two; + if ((Y - One == Zero)) StickyBit = S; + } + else StickyBit = S; + } + } + } + } + } + } + if (StickyBit == One) printf("Sticky bit apparently used correctly.\n"); + else printf("Sticky bit used incorrectly or not at all.\n"); + TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No || + RMult == Other || RDiv == Other || RAddSub == Other), + "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"); + /*=============================================*/ + Milestone = 60; + /*=============================================*/ + printf("\n"); + printf("Does Multiplication commute? "); + printf("Testing on %d random pairs.\n", NoTrials); + Random9 = SQRT(3.0); + Random1 = Third; + I = 1; + do { + X = Random(); + Y = Random(); + Z9 = Y * X; + Z = X * Y; + Z9 = Z - Z9; + I = I + 1; + } while ( ! ((I > NoTrials) || (Z9 != Zero))); + if (I == NoTrials) { + Random1 = One + Half / Three; + Random2 = (U2 + U1) + One; + Z = Random1 * Random2; + Y = Random2 * Random1; + Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half / + Three) * ((U2 + U1) + One); + } + if (! ((I == NoTrials) || (Z9 == Zero))) + BadCond(Defect, "X * Y == Y * X trial fails.\n"); + else printf(" No failures found in %d integer pairs.\n", NoTrials); + /*=============================================*/ + Milestone = 70; + /*=============================================*/ + printf("\nRunning test of square root(x).\n"); + TstCond (Failure, (Zero == SQRT(Zero)) + && (- Zero == SQRT(- Zero)) + && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong"); + MinSqEr = Zero; + MaxSqEr = Zero; + J = Zero; + X = Radix; + OneUlp = U2; + SqXMinX (Serious); + X = BInvrse; + OneUlp = BInvrse * U1; + SqXMinX (Serious); + X = U1; + OneUlp = U1 * U1; + SqXMinX (Serious); + if (J != Zero) Pause(); + printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials); + J = Zero; + X = Two; + Y = Radix; + if ((Radix != One)) do { + X = Y; + Y = Radix * Y; + } while ( ! ((Y - X >= NoTrials))); + OneUlp = X * U2; + I = 1; + while (I <= NoTrials) { + X = X + One; + SqXMinX (Defect); + if (J > Zero) break; + I = I + 1; + } + printf("Test for sqrt monotonicity.\n"); + I = - 1; + X = BMinusU2; + Y = Radix; + Z = Radix + Radix * U2; + NotMonot = False; + Monot = False; + while ( ! (NotMonot || Monot)) { + I = I + 1; + X = SQRT(X); + Q = SQRT(Y); + Z = SQRT(Z); + if ((X > Q) || (Q > Z)) NotMonot = True; + else { + Q = FLOOR(Q + Half); + if ((I > 0) || (Radix == Q * Q)) Monot = True; + else if (I > 0) { + if (I > 1) Monot = True; + else { + Y = Y * BInvrse; + X = Y - U1; + Z = Y + U1; + } + } + else { + Y = Q; + X = Y - U2; + Z = Y + U2; + } + } + } + if (Monot) printf("sqrt has passed a test for Monotonicity.\n"); + else { + BadCond(Defect, ""); + printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y); + } + /*=============================================*/ + /*SPLIT + } +#include "paranoia.h" +part5(){ +*/ + Milestone = 80; + /*=============================================*/ + MinSqEr = MinSqEr + Half; + MaxSqEr = MaxSqEr - Half; + Y = (SQRT(One + U2) - One) / U2; + SqEr = (Y - One) + U2 / Eight; + if (SqEr > MaxSqEr) MaxSqEr = SqEr; + SqEr = Y + U2 / Eight; + if (SqEr < MinSqEr) MinSqEr = SqEr; + Y = ((SQRT(F9) - U2) - (One - U2)) / U1; + SqEr = Y + U1 / Eight; + if (SqEr > MaxSqEr) MaxSqEr = SqEr; + SqEr = (Y + One) + U1 / Eight; + if (SqEr < MinSqEr) MinSqEr = SqEr; + OneUlp = U2; + X = OneUlp; + for( Indx = 1; Indx <= 3; ++Indx) { + Y = SQRT((X + U1 + X) + F9); + Y = ((Y - U2) - ((One - U2) + X)) / OneUlp; + Z = ((U1 - X) + F9) * Half * X * X / OneUlp; + SqEr = (Y + Half) + Z; + if (SqEr < MinSqEr) MinSqEr = SqEr; + SqEr = (Y - Half) + Z; + if (SqEr > MaxSqEr) MaxSqEr = SqEr; + if (((Indx == 1) || (Indx == 3))) + X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp))); + else { + OneUlp = U1; + X = - OneUlp; + } + } + /*=============================================*/ + Milestone = 85; + /*=============================================*/ + SqRWrng = False; + Anomaly = False; + RSqrt = Other; /* ~dgh */ + if (Radix != One) { + printf("Testing whether sqrt is rounded or chopped.\n"); + D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision))); + /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */ + X = D / Radix; + Y = D / A1; + if ((X != FLOOR(X)) || (Y != FLOOR(Y))) { + Anomaly = True; + } + else { + X = Zero; + Z2 = X; + Y = One; + Y2 = Y; + Z1 = Radix - One; + FourD = Four * D; + do { + if (Y2 > Z2) { + Q = Radix; + Y1 = Y; + do { + X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1); + Q = Y1; + Y1 = X1; + } while ( ! (X1 <= Zero)); + if (Q <= One) { + Z2 = Y2; + Z = Y; + } + } + Y = Y + Two; + X = X + Eight; + Y2 = Y2 + X; + if (Y2 >= FourD) Y2 = Y2 - FourD; + } while ( ! (Y >= D)); + X8 = FourD - Z2; + Q = (X8 + Z * Z) / FourD; + X8 = X8 / Eight; + if (Q != FLOOR(Q)) Anomaly = True; + else { + Break = False; + do { + X = Z1 * Z; + X = X - FLOOR(X / Radix) * Radix; + if (X == One) + Break = True; + else + Z1 = Z1 - One; + } while ( ! (Break || (Z1 <= Zero))); + if ((Z1 <= Zero) && (! Break)) Anomaly = True; + else { + if (Z1 > RadixD2) Z1 = Z1 - Radix; + do { + NewD(); + } while ( ! (U2 * D >= F9)); + if (D * Radix - D != W - D) Anomaly = True; + else { + Z2 = D; + I = 0; + Y = D + (One + Z) * Half; + X = D + Z + Q; + SR3750(); + Y = D + (One - Z) * Half + D; + X = D - Z + D; + X = X + Q + X; + SR3750(); + NewD(); + if (D - Z2 != W - Z2) Anomaly = True; + else { + Y = (D - Z2) + (Z2 + (One - Z) * Half); + X = (D - Z2) + (Z2 - Z + Q); + SR3750(); + Y = (One + Z) * Half; + X = Q; + SR3750(); + if (I == 0) Anomaly = True; + } + } + } + } + } + if ((I == 0) || Anomaly) { + BadCond(Failure, "Anomalous arithmetic with Integer < "); + printf("Radix^Precision = %.7e\n", W); + printf(" fails test whether sqrt rounds or chops.\n"); + SqRWrng = True; + } + } + if (! Anomaly) { + if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) { + RSqrt = Rounded; + printf("Square root appears to be correctly rounded.\n"); + } + else { + if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half) + || (MinSqEr + Radix < Half)) SqRWrng = True; + else { + RSqrt = Chopped; + printf("Square root appears to be chopped.\n"); + } + } + } + if (SqRWrng) { + printf("Square root is neither chopped nor correctly rounded.\n"); + printf("Observed errors run from %.7e ", MinSqEr - Half); + printf("to %.7e ulps.\n", Half + MaxSqEr); + TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix, + "sqrt gets too many last digits wrong"); + } + /*=============================================*/ + Milestone = 90; + /*=============================================*/ + Pause(); + printf("Testing powers Z^i for small Integers Z and i.\n"); + N = 0; + /* ... test powers of zero. */ + I = 0; + Z = -Zero; + M = 3.0; + Break = False; + do { + X = One; + SR3980(); + if (I <= 10) { + I = 1023; + SR3980(); + } + if (Z == MinusOne) Break = True; + else { + Z = MinusOne; + PrintIfNPositive(); + N = 0; + /* .. if(-1)^N is invalid, replace MinusOne by One. */ + I = - 4; + } + } while ( ! Break); + PrintIfNPositive(); + N1 = N; + N = 0; + Z = A1; + M = FLOOR(Two * LOG(W) / LOG(A1)); + Break = False; + do { + X = Z; + I = 1; + SR3980(); + if (Z == AInvrse) Break = True; + else Z = AInvrse; + } while ( ! (Break)); + /*=============================================*/ + Milestone = 100; + /*=============================================*/ + /* Powers of Radix have been tested, */ + /* next try a few primes */ + M = NoTrials; + Z = Three; + do { + X = Z; + I = 1; + SR3980(); + do { + Z = Z + Two; + } while ( Three * FLOOR(Z / Three) == Z ); + } while ( Z < Eight * Three ); + if (N > 0) { + printf("Errors like this may invalidate financial calculations\n"); + printf("\tinvolving interest rates.\n"); + } + PrintIfNPositive(); + N += N1; + if (N == 0) printf("... no discrepancis found.\n"); + if (N > 0) Pause(); + else printf("\n"); + /*=============================================*/ + /*SPLIT + } +#include "paranoia.h" +part6(){ +*/ + Milestone = 110; + /*=============================================*/ + printf("Seeking Underflow thresholds UfThold and E0.\n"); + D = U1; + if (Precision != FLOOR(Precision)) { + D = BInvrse; + X = Precision; + do { + D = D * BInvrse; + X = X - One; + } while ( X > Zero); + } + Y = One; + Z = D; + /* ... D is power of 1/Radix < 1. */ + do { + C = Y; + Y = Z; + Z = Y * Y; + } while ((Y > Z) && (Z + Z > Z)); + Y = C; + Z = Y * D; + do { + C = Y; + Y = Z; + Z = Y * D; + } while ((Y > Z) && (Z + Z > Z)); + if (Radix < Two) HInvrse = Two; + else HInvrse = Radix; + H = One / HInvrse; + /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */ + CInvrse = One / C; + E0 = C; + Z = E0 * H; + /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */ + do { + Y = E0; + E0 = Z; + Z = E0 * H; + } while ((E0 > Z) && (Z + Z > Z)); + UfThold = E0; + E1 = Zero; + Q = Zero; + E9 = U2; + S = One + E9; + D = C * S; + if (D <= C) { + E9 = Radix * U2; + S = One + E9; + D = C * S; + if (D <= C) { + BadCond(Failure, "multiplication gets too many last digits wrong.\n"); + Underflow = E0; + Y1 = Zero; + PseudoZero = Z; + Pause(); + } + } + else { + Underflow = D; + PseudoZero = Underflow * H; + UfThold = Zero; + do { + Y1 = Underflow; + Underflow = PseudoZero; + if (E1 + E1 <= E1) { + Y2 = Underflow * HInvrse; + E1 = FABS(Y1 - Y2); + Q = Y1; + if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1; + } + PseudoZero = PseudoZero * H; + } while ((Underflow > PseudoZero) + && (PseudoZero + PseudoZero > PseudoZero)); + } + /* Comment line 4530 .. 4560 */ + if (PseudoZero != Zero) { + printf("\n"); + Z = PseudoZero; + /* ... Test PseudoZero for "phoney- zero" violates */ + /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero + ... */ + if (PseudoZero <= Zero) { + BadCond(Failure, "Positive expressions can underflow to an\n"); + printf("allegedly negative value\n"); + printf("PseudoZero that prints out as: %g .\n", PseudoZero); + X = - PseudoZero; + if (X <= Zero) { + printf("But -PseudoZero, which should be\n"); + printf("positive, isn't; it prints out as %g .\n", X); + } + } + else { + BadCond(Flaw, "Underflow can stick at an allegedly positive\n"); + printf("value PseudoZero that prints out as %g .\n", PseudoZero); + } + TstPtUf(); + } + /*=============================================*/ + Milestone = 120; + /*=============================================*/ + if (CInvrse * Y > CInvrse * Y1) { + S = H * S; + E0 = Underflow; + } + if (! ((E1 == Zero) || (E1 == E0))) { + BadCond(Defect, ""); + if (E1 < E0) { + printf("Products underflow at a higher"); + printf(" threshold than differences.\n"); + if (PseudoZero == Zero) + E0 = E1; + } + else { + printf("Difference underflows at a higher"); + printf(" threshold than products.\n"); + } + } + printf("Smallest strictly positive number found is E0 = %g .\n", E0); + Z = E0; + TstPtUf(); + Underflow = E0; + if (N == 1) Underflow = Y; + I = 4; + if (E1 == Zero) I = 3; + if (UfThold == Zero) I = I - 2; + UfNGrad = True; + switch (I) { + case 1: + UfThold = Underflow; + if ((CInvrse * Q) != ((CInvrse * Y) * S)) { + UfThold = Y; + BadCond(Failure, "Either accuracy deteriorates as numbers\n"); + printf("approach a threshold = %.17e\n", UfThold);; + printf(" coming down from %.17e\n", C); + printf(" or else multiplication gets too many last digits wrong.\n"); + } + Pause(); + break; + + case 2: + BadCond(Failure, "Underflow confuses Comparison, which alleges that\n"); + printf("Q == Y while denying that |Q - Y| == 0; these values\n"); + printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2); + printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2)); + UfThold = Q; + break; + + case 3: + X = X; + break; + + case 4: + if ((Q == UfThold) && (E1 == E0) + && (FABS( UfThold - E1 / E9) <= E1)) { + UfNGrad = False; + printf("Underflow is gradual; it incurs Absolute Error =\n"); + printf("(roundoff in UfThold) < E0.\n"); + Y = E0 * CInvrse; + Y = Y * (OneAndHalf + U2); + X = CInvrse * (One + U2); + Y = Y / X; + IEEE = (Y == E0); + } + } + if (UfNGrad) { + printf("\n"); + sigsave = sigfpe; + if (setjmp(ovfl_buf)) { + printf("Underflow / UfThold failed!\n"); + R = H + H; + } + else R = SQRT(Underflow / UfThold); + sigsave = 0; + if (R <= H) { + Z = R * UfThold; + X = Z * (One + R * H * (One + H)); + } + else { + Z = UfThold; + X = Z * (One + H * H * (One + H)); + } + if (! ((X == Z) || (X - Z != Zero))) { + BadCond(Flaw, ""); + printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z); + Z9 = X - Z; + printf("yet X - Z yields %.17e .\n", Z9); + printf(" Should this NOT signal Underflow, "); + printf("this is a SERIOUS DEFECT\nthat causes "); + printf("confusion when innocent statements like\n");; + printf(" if (X == Z) ... else"); + printf(" ... (f(X) - f(Z)) / (X - Z) ...\n"); + printf("encounter Division by Zero although actually\n"); + sigsave = sigfpe; + if (setjmp(ovfl_buf)) printf("X / Z fails!\n"); + else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half); + sigsave = 0; + } + } + printf("The Underflow threshold is %.17e, %s\n", UfThold, + " below which"); + printf("calculation may suffer larger Relative error than "); + printf("merely roundoff.\n"); + Y2 = U1 * U1; + Y = Y2 * Y2; + Y2 = Y * U1; + if (Y2 <= UfThold) { + if (Y > E0) { + BadCond(Defect, ""); + I = 5; + } + else { + BadCond(Serious, ""); + I = 4; + } + printf("Range is too narrow; U1^%d Underflows.\n", I); + } + /*=============================================*/ + /*SPLIT + } +#include "paranoia.h" +part7(){ +*/ + Milestone = 130; + /*=============================================*/ + Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty; + Y2 = Y + Y; + printf("Since underflow occurs below the threshold\n"); + printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y); + printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y); + V9 = POW(HInvrse, Y2); + printf("actually calculating yields: %.17e .\n", V9); + if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) { + BadCond(Serious, "this is not between 0 and underflow\n"); + printf(" threshold = %.17e .\n", UfThold); + } + else if (! (V9 > UfThold * (One + E9))) + printf("This computed value is O.K.\n"); + else { + BadCond(Defect, "this is not between 0 and underflow\n"); + printf(" threshold = %.17e .\n", UfThold); + } + /*=============================================*/ + Milestone = 140; + /*=============================================*/ + printf("\n"); + /* ...calculate Exp2 == exp(2) == 7.389056099... */ + X = Zero; + I = 2; + Y = Two * Three; + Q = Zero; + N = 0; + do { + Z = X; + I = I + 1; + Y = Y / (I + I); + R = Y + Q; + X = Z + R; + Q = (Z - X) + R; + } while(X > Z); + Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo); + X = Z * Z; + Exp2 = X * X; + X = F9; + Y = X - U1; + printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n", + Exp2); + for(I = 1;;) { + Z = X - BInvrse; + Z = (X + One) / (Z - (One - BInvrse)); + Q = POW(X, Z) - Exp2; + if (FABS(Q) > TwoForty * U2) { + N = 1; + V9 = (X - BInvrse) - (One - BInvrse); + BadCond(Defect, "Calculated"); + printf(" %.17e for\n", POW(X,Z)); + printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z); + printf("\tdiffers from correct value by %.17e .\n", Q); + printf("\tThis much error may spoil financial\n"); + printf("\tcalculations involving tiny interest rates.\n"); + break; + } + else { + Z = (Y - X) * Two + Y; + X = Y; + Y = Z; + Z = One + (X - F9)*(X - F9); + if (Z > One && I < NoTrials) I++; + else { + if (X > One) { + if (N == 0) + printf("Accuracy seems adequate.\n"); + break; + } + else { + X = One + U2; + Y = U2 + U2; + Y += X; + I = 1; + } + } + } + } + /*=============================================*/ + Milestone = 150; + /*=============================================*/ + printf("Testing powers Z^Q at four nearly extreme values.\n"); + N = 0; + Z = A1; + Q = FLOOR(Half - LOG(C) / LOG(A1)); + Break = False; + do { + X = CInvrse; + Y = POW(Z, Q); + IsYeqX(); + Q = - Q; + X = C; + Y = POW(Z, Q); + IsYeqX(); + if (Z < One) Break = True; + else Z = AInvrse; + } while ( ! (Break)); + PrintIfNPositive(); + if (N == 0) printf(" ... no discrepancies found.\n"); + printf("\n"); + + /*=============================================*/ + Milestone = 160; + /*=============================================*/ + Pause(); + printf("Searching for Overflow threshold:\n"); + printf("This may generate an error.\n"); + Y = - CInvrse; + V9 = HInvrse * Y; + sigsave = sigfpe; + if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; } + do { + V = Y; + Y = V9; + V9 = HInvrse * Y; + } while(V9 < Y); + I = 1; +overflow: + sigsave = 0; + Z = V9; + printf("Can `Z = -Y' overflow?\n"); + printf("Trying it on Y = %.17e .\n", Y); + V9 = - Y; + V0 = V9; + if (V - Y == V + V0) printf("Seems O.K.\n"); + else { + printf("finds a "); + BadCond(Flaw, "-(-Y) differs from Y.\n"); + } + if (Z != Y) { + BadCond(Serious, ""); + printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z); + } + if (I) { + Y = V * (HInvrse * U2 - HInvrse); + Z = Y + ((One - HInvrse) * U2) * V; + if (Z < V0) Y = Z; + if (Y < V0) V = Y; + if (V0 - V < V0) V = V0; + } + else { + V = Y * (HInvrse * U2 - HInvrse); + V = V + ((One - HInvrse) * U2) * Y; + } + printf("Overflow threshold is V = %.17e .\n", V); + if (I) printf("Overflow saturates at V0 = %.17e .\n", V0); + else printf("There is no saturation value because the system traps on overflow.\n"); + V9 = V * One; + printf("No Overflow should be signaled for V * 1 = %.17e\n", V9); + V9 = V / One; + printf(" nor for V / 1 = %.17e .\n", V9); + printf("Any overflow signal separating this * from the one\n"); + printf("above is a DEFECT.\n"); + /*=============================================*/ + Milestone = 170; + /*=============================================*/ + if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) { + BadCond(Failure, "Comparisons involving "); + printf("+-%g, +-%g\nand +-%g are confused by Overflow.", + V, V0, UfThold); + } + /*=============================================*/ + Milestone = 175; + /*=============================================*/ + printf("\n"); + for(Indx = 1; Indx <= 3; ++Indx) { + switch (Indx) { + case 1: Z = UfThold; break; + case 2: Z = E0; break; + case 3: Z = PseudoZero; break; + } + if (Z != Zero) { + V9 = SQRT(Z); + Y = V9 * V9; + if (Y / (One - Radix * E9) < Z + || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */ + if (V9 > U1) BadCond(Serious, ""); + else BadCond(Defect, ""); + printf("Comparison alleges that what prints as Z = %.17e\n", Z); + printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y); + } + } + } + /*=============================================*/ + Milestone = 180; + /*=============================================*/ + for(Indx = 1; Indx <= 2; ++Indx) { + if (Indx == 1) Z = V; + else Z = V0; + V9 = SQRT(Z); + X = (One - Radix * E9) * V9; + V9 = V9 * X; + if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) { + Y = V9; + if (X < W) BadCond(Serious, ""); + else BadCond(Defect, ""); + printf("Comparison alleges that Z = %17e\n", Z); + printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y); + } + } + /*=============================================*/ + /*SPLIT + } +#include "paranoia.h" +part8(){ +*/ + Milestone = 190; + /*=============================================*/ + Pause(); + X = UfThold * V; + Y = Radix * Radix; + if (X*Y < One || X > Y) { + if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly"); + else BadCond(Flaw, ""); + + printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n", + X, "is too far from 1.\n"); + } + /*=============================================*/ + Milestone = 200; + /*=============================================*/ + for (Indx = 1; Indx <= 5; ++Indx) { + X = F9; + switch (Indx) { + case 2: X = One + U2; break; + case 3: X = V; break; + case 4: X = UfThold; break; + case 5: X = Radix; + } + Y = X; + sigsave = sigfpe; + if (setjmp(ovfl_buf)) + printf(" X / X traps when X = %g\n", X); + else { + V9 = (Y / X - Half) - Half; + if (V9 == Zero) continue; + if (V9 == - U1 && Indx < 5) BadCond(Flaw, ""); + else BadCond(Serious, ""); + printf(" X / X differs from 1 when X = %.17e\n", X); + printf(" instead, X / X - 1/2 - 1/2 = %.17e .\n", V9); + } + sigsave = 0; + } + /*=============================================*/ + Milestone = 210; + /*=============================================*/ + MyZero = Zero; + printf("\n"); + printf("What message and/or values does Division by Zero produce?\n") ; +#ifndef NOPAUSE + printf("This can interupt your program. You can "); + printf("skip this part if you wish.\n"); + printf("Do you wish to compute 1 / 0? "); + fflush(stdout); + read (KEYBOARD, ch, 8); + if ((ch[0] == 'Y') || (ch[0] == 'y')) { +#endif + sigsave = sigfpe; + printf(" Trying to compute 1 / 0 produces ..."); + if (!setjmp(ovfl_buf)) printf(" %.7e .\n", One / MyZero); + sigsave = 0; +#ifndef NOPAUSE + } + else printf("O.K.\n"); + printf("\nDo you wish to compute 0 / 0? "); + fflush(stdout); + read (KEYBOARD, ch, 80); + if ((ch[0] == 'Y') || (ch[0] == 'y')) { +#endif + sigsave = sigfpe; + printf("\n Trying to compute 0 / 0 produces ..."); + if (!setjmp(ovfl_buf)) printf(" %.7e .\n", Zero / MyZero); + sigsave = 0; +#ifndef NOPAUSE + } + else printf("O.K.\n"); +#endif + /*=============================================*/ + Milestone = 220; + /*=============================================*/ + Pause(); + printf("\n"); + { + static char *msg[] = { + "FAILUREs encountered =", + "SERIOUS DEFECTs discovered =", + "DEFECTs discovered =", + "FLAWs discovered =" }; + int i; + for(i = 0; i < 4; i++) if (ErrCnt[i]) + printf("The number of %-29s %d.\n", + msg[i], ErrCnt[i]); + } + printf("\n"); + if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect] + + ErrCnt[Flaw]) > 0) { + if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[ + Defect] == 0) && (ErrCnt[Flaw] > 0)) { + printf("The arithmetic diagnosed seems "); + printf("Satisfactory though flawed.\n"); + } + if ((ErrCnt[Failure] + ErrCnt[Serious] == 0) + && ( ErrCnt[Defect] > 0)) { + printf("The arithmetic diagnosed may be Acceptable\n"); + printf("despite inconvenient Defects.\n"); + } + if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) { + printf("The arithmetic diagnosed has "); + printf("unacceptable Serious Defects.\n"); + } + if (ErrCnt[Failure] > 0) { + printf("Potentially fatal FAILURE may have spoiled this"); + printf(" program's subsequent diagnoses.\n"); + } + } + else { + printf("No failures, defects nor flaws have been discovered.\n"); + if (! ((RMult == Rounded) && (RDiv == Rounded) + && (RAddSub == Rounded) && (RSqrt == Rounded))) + printf("The arithmetic diagnosed seems Satisfactory.\n"); + else { + if (StickyBit >= One && + (Radix - Two) * (Radix - Nine - One) == Zero) { + printf("Rounding appears to conform to "); + printf("the proposed IEEE standard P"); + if ((Radix == Two) && + ((Precision - Four * Three * Two) * + ( Precision - TwentySeven - + TwentySeven + One) == Zero)) + printf("754"); + else printf("854"); + if (IEEE) printf(".\n"); + else { + printf(",\nexcept for possibly Double Rounding"); + printf(" during Gradual Underflow.\n"); + } + } + printf("The arithmetic diagnosed appears to be Excellent!\n"); + } + } + if (fpecount) + printf("\nA total of %d floating point exceptions were registered.\n", + fpecount); + printf("END OF TEST.\n"); + return 0; + } + +/*SPLIT subs.c +#include "paranoia.h" +*/ + +/* Sign */ + +FLOAT Sign (X) +FLOAT X; +{ return X >= 0. ? 1.0 : -1.0; } + +/* Pause */ + +Pause() +{ +#ifndef NOPAUSE + char ch[8]; + + printf("\nTo continue, press RETURN"); + fflush(stdout); + read(KEYBOARD, ch, 8); +#endif + printf("\nDiagnosis resumes after milestone Number %d", Milestone); + printf(" Page: %d\n\n", PageNo); + ++Milestone; + ++PageNo; + } + + /* TstCond */ + +TstCond (K, Valid, T) +int K, Valid; +char *T; +{ if (! Valid) { BadCond(K,T); printf(".\n"); } } + +BadCond(K, T) +int K; +char *T; +{ + static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" }; + + ErrCnt [K] = ErrCnt [K] + 1; + printf("%s: %s", msg[K], T); + } + +/* Random */ +/* Random computes + X = (Random1 + Random9)^5 + Random1 = X - FLOOR(X) + 0.000005 * X; + and returns the new value of Random1 +*/ + +FLOAT Random() +{ + FLOAT X, Y; + + X = Random1 + Random9; + Y = X * X; + Y = Y * Y; + X = X * Y; + Y = X - FLOOR(X); + Random1 = Y + X * 0.000005; + return(Random1); + } + +/* SqXMinX */ + +SqXMinX (ErrKind) +int ErrKind; +{ + FLOAT XA, XB; + + XB = X * BInvrse; + XA = X - XB; + SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp; + if (SqEr != Zero) { + if (SqEr < MinSqEr) MinSqEr = SqEr; + if (SqEr > MaxSqEr) MaxSqEr = SqEr; + J = J + 1.0; + BadCond(ErrKind, "\n"); + printf("sqrt( %.17e) - %.17e = %.17e\n", X * X, X, OneUlp * SqEr); + printf("\tinstead of correct value 0 .\n"); + } + } + +/* NewD */ + +NewD() +{ + X = Z1 * Q; + X = FLOOR(Half - X / Radix) * Radix + X; + Q = (Q - X * Z) / Radix + X * X * (D / Radix); + Z = Z - Two * X * D; + if (Z <= Zero) { + Z = - Z; + Z1 = - Z1; + } + D = Radix * D; + } + +/* SR3750 */ + +SR3750() +{ + if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) { + I = I + 1; + X2 = SQRT(X * D); + Y2 = (X2 - Z2) - (Y - Z2); + X2 = X8 / (Y - Half); + X2 = X2 - Half * X2 * X2; + SqEr = (Y2 + Half) + (Half - X2); + if (SqEr < MinSqEr) MinSqEr = SqEr; + SqEr = Y2 - X2; + if (SqEr > MaxSqEr) MaxSqEr = SqEr; + } + } + +/* IsYeqX */ + +IsYeqX() +{ + if (Y != X) { + if (N <= 0) { + if (Z == Zero && Q <= Zero) + printf("WARNING: computing\n"); + else BadCond(Defect, "computing\n"); + printf("\t(%.17e) ^ (%.17e)\n", Z, Q); + printf("\tyielded %.17e;\n", Y); + printf("\twhich compared unequal to correct %.17e ;\n", + X); + printf("\t\tthey differ by %.17e .\n", Y - X); + } + N = N + 1; /* ... count discrepancies. */ + } + } + +/* SR3980 */ + +SR3980() +{ + do { + Q = (FLOAT) I; + Y = POW(Z, Q); + IsYeqX(); + if (++I > M) break; + X = Z * X; + } while ( X < W ); + } + +/* PrintIfNPositive */ + +PrintIfNPositive() +{ + if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N); + } + +/* TstPtUf */ + +TstPtUf() +{ + N = 0; + if (Z != Zero) { + printf("Since comparison denies Z = 0, evaluating "); + printf("(Z + Z) / Z should be safe.\n"); + sigsave = sigfpe; + if (setjmp(ovfl_buf)) goto very_serious; + Q9 = (Z + Z) / Z; + printf("What the machine gets for (Z + Z) / Z is %.17e .\n", + Q9); + if (FABS(Q9 - Two) < Radix * U2) { + printf("This is O.K., provided Over/Underflow"); + printf(" has NOT just been signaled.\n"); + } + else { + if ((Q9 < One) || (Q9 > Two)) { +very_serious: + N = 1; + ErrCnt [Serious] = ErrCnt [Serious] + 1; + printf("This is a VERY SERIOUS DEFECT!\n"); + } + else { + N = 1; + ErrCnt [Defect] = ErrCnt [Defect] + 1; + printf("This is a DEFECT!\n"); + } + } + sigsave = 0; + V9 = Z * One; + Random1 = V9; + V9 = One * Z; + Random2 = V9; + V9 = Z / One; + if ((Z == Random1) && (Z == Random2) && (Z == V9)) { + if (N > 0) Pause(); + } + else { + N = 1; + BadCond(Defect, "What prints as Z = "); + printf("%.17e\n\tcompares different from ", Z); + if (Z != Random1) printf("Z * 1 = %.17e ", Random1); + if (! ((Z == Random2) + || (Random2 == Random1))) + printf("1 * Z == %g\n", Random2); + if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9); + if (Random2 != Random1) { + ErrCnt [Defect] = ErrCnt [Defect] + 1; + BadCond(Defect, "Multiplication does not commute!\n"); + printf("\tComparison alleges that 1 * Z = %.17e\n", + Random2); + printf("\tdiffers from Z * 1 = %.17e\n", Random1); + } + Pause(); + } + } + } + +notify(s) +char *s; +{ + printf("%s test appears to be inconsistent...\n", s); + printf(" PLEASE NOTIFY KARPINKSI!\n"); + } + +/*SPLIT msgs.c */ + +/* Instructions */ + +msglist(s) +char **s; +{ while(*s) printf("%s\n", *s++); } + +Instructions() +{ + static char *instr[] = { + "Lest this program stop prematurely, i.e. before displaying\n", + " `END OF TEST',\n", + "try to persuade the computer NOT to terminate execution when an", + "error like Over/Underflow or Division by Zero occurs, but rather", + "to persevere with a surrogate value after, perhaps, displaying some", + "warning. If persuasion avails naught, don't despair but run this", + "program anyway to see how many milestones it passes, and then", + "amend it to make further progress.\n", + "Answer questions with Y, y, N or n (unless otherwise indicated).\n", + 0}; + + msglist(instr); + } + +/* Heading */ + +Heading() +{ + static char *head[] = { + "Users are invited to help debug and augment this program so it will", + "cope with unanticipated and newly uncovered arithmetic pathologies.\n", + "Please send suggestions and interesting results to", + "\tRichard Karpinski", + "\tComputer Center U-76", + "\tUniversity of California", + "\tSan Francisco, CA 94143-0704, USA\n", + "In doing so, please include the following information:", +#ifdef Single + "\tPrecision:\tsingle;", +#else + "\tPrecision:\tdouble;", +#endif + "\tVersion:\t10 February 1989;", + "\tComputer:\n", + "\tCompiler:\n", + "\tOptimization level:\n", + "\tOther relevant compiler options:", + 0}; + + msglist(head); + } + +/* Characteristics */ + +Characteristics() +{ + static char *chars[] = { + "Running this program should reveal these characteristics:", + " Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...", + " Precision = number of significant digits carried.", + " U2 = Radix/Radix^Precision = One Ulp", + "\t(OneUlpnit in the Last Place) of 1.000xxx .", + " U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .", + " Adequacy of guard digits for Mult., Div. and Subt.", + " Whether arithmetic is chopped, correctly rounded, or something else", + "\tfor Mult., Div., Add/Subt. and Sqrt.", + " Whether a Sticky Bit used correctly for rounding.", + " UnderflowThreshold = an underflow threshold.", + " E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.", + " V = an overflow threshold, roughly.", + " V0 tells, roughly, whether Infinity is represented.", + " Comparisions are checked for consistency with subtraction", + "\tand for contamination with pseudo-zeros.", + " Sqrt is tested. Y^X is not tested.", + " Extra-precise subexpressions are revealed but NOT YET tested.", + " Decimal-Binary conversion is NOT YET tested for accuracy.", + 0}; + + msglist(chars); + } + +History() + +{ /* History */ + /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner, + with further massaging by David M. Gay. */ + + static char *hist[] = { + "The program attempts to discriminate among", + " FLAWs, like lack of a sticky bit,", + " Serious DEFECTs, like lack of a guard digit, and", + " FAILUREs, like 2+2 == 5 .", + "Failures may confound subsequent diagnoses.\n", + "The diagnostic capabilities of this program go beyond an earlier", + "program called `MACHAR', which can be found at the end of the", + "book `Software Manual for the Elementary Functions' (1980) by", + "W. J. Cody and W. Waite. Although both programs try to discover", + "the Radix, Precision and range (over/underflow thresholds)", + "of the arithmetic, this program tries to cope with a wider variety", + "of pathologies, and to say how well the arithmetic is implemented.", + "\nThe program is based upon a conventional radix representation for", + "floating-point numbers, but also allows logarithmic encoding", + "as used by certain early WANG machines.\n", + "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;", + "see source comments for more history.", + 0}; + + msglist(hist); + } + +double +pow(x, y) /* return x ^ y (exponentiation) */ +double x, y; +{ + extern double exp(), frexp(), ldexp(), log(), modf(); + double xy, ye; + long i; + int ex, ey = 0, flip = 0; + + if (!y) return 1.0; + + if ((y < -1100. || y > 1100.) && x != -1.) return exp(y * log(x)); + + if (y < 0.) { y = -y; flip = 1; } + y = modf(y, &ye); + if (y) xy = exp(y * log(x)); + else xy = 1.0; + /* next several lines assume >= 32 bit integers */ + x = frexp(x, &ex); + if (i = ye) for(;;) { + if (i & 1) { xy *= x; ey += ex; } + if (!(i >>= 1)) break; + x *= x; + ex *= 2; + if (x < .5) { x *= 2.; ex -= 1; } + } + if (flip) { xy = 1. / xy; ey = -ey; } + return ldexp(xy, ey); +} |