diff options
Diffstat (limited to 'libm/ldouble/lparanoi.c')
| -rw-r--r-- | libm/ldouble/lparanoi.c | 2348 |
1 files changed, 2348 insertions, 0 deletions
diff --git a/libm/ldouble/lparanoi.c b/libm/ldouble/lparanoi.c new file mode 100644 index 000000000..eb8fd25c7 --- /dev/null +++ b/libm/ldouble/lparanoi.c @@ -0,0 +1,2348 @@ +/* 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 Mountainn 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). + +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! = .*!;! +/^Guard/,/^Round/s/^/extern / +/^jmp_buf/s/^/extern / +/^Sig_type/s/^/extern / +a +extern int sigfpe(); +. +w paranoia.h +q + +*/ + +#include <stdio.h> +#ifndef NOSIGNAL +#include <signal.h> +#endif +#include <setjmp.h> + +#define Ldouble +/*#define Single*/ + +#ifdef Single +#define NPRT 2 +extern double fabs(), floor(), log(), pow(), sqrt(); +#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)) +#define FSETUP sprec +/*sprec() { }*/ +#else +#ifdef Ldouble +#define NPRT 6 +extern long double fabsl(), floorl(), logl(), powl(), sqrtl(); +#define FLOAT long double +#define FABS(x) fabsl(x) +#define FLOOR(x) floorl(x) +#define LOG(x) logl(x) +#define POW(x,y) powl(x,y) +#define SQRT(x) sqrtl(x) +#define FSETUP ldprec +#else +#define NPRT 4 +extern double fabs(), floor(), log(), pow(), sqrt(); +#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) +/*double __sqrtdf2(); +#define SQRT(x) __sqrtdf2(x) +*/ +#define FSETUP dprec +/* dprec() { } */ +#endif +#endif + +jmp_buf ovfl_buf; +typedef int (*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; +static 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; +static 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 */ +sigfpe() +{ + fpecount++; + printf("\n* * * FLOATING-POINT ERROR * * *\n"); + fflush(stdout); + if (sigsave) { +#ifndef NOSIGNAL + signal(SIGFPE, sigsave); +#endif + sigsave = 0; + longjmp(ovfl_buf, 1); + } + abort(); +} + + +FLOAT Ptemp; + +pnum( x ) +FLOAT *x; +{ +char str[30]; +double d; +unsigned short *p; +int i; + +p = (unsigned short *)x; +for( i=0; i<NPRT; i++ ) + printf( "%04x ", *p++ & 0xffff ); +#ifdef Ldouble +e64toasc( x, str, 20 ); +#else +#ifdef Single +e24toasc( x, str, 20 ); +#else +e53toasc( x, str, 20 ); +#endif +#endif +printf( " = %s\n", str ); +/* +d = *x; +printf( " = %.16e\n", d ); +*/ +} + + + +main() +{ +/* noexcept(); */ + FSETUP(); + /* 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) { + U1 = 0.001; + Radix = 1; + TstPtUf(); + } + else { + ErrCnt[Failure] = ErrCnt[Failure] + 1; + printf("Comparison alleges that -0.0 is Non-zero!\n"); + } + 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 = " ); + pnum( &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 = " ); + pnum( &U1 ); + printf("U2 = "); + pnum( &U2 ); + printf("Recalculating radix and precision."); + + /*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 = " ); + pnum( &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 = " ); + pnum( &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 " ); + pnum( &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 " ); + pnum( &X1 ); + pnum( &Y1 ); + pnum( &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 = " ); + pnum( &U1 ); + printf( "Z1 - U1 = " ); + Ptemp = Z1-U1; + pnum( &Ptemp ); + printf("\tU2 = " ); + pnum( &U2 ); + Ptemp = Z2-U2; + printf( "Z2 - U2 = " ); + pnum( &Ptemp ); + } + else { + if ((Z1 <= Zero) || (Z2 <= Zero)) { + printf("Because of unusual Radix = "); + pnum( &Radix ); + printf(", or exact rational arithmetic a result\n"); + printf("Z1 = " ); + pnum( &Z1 ); + printf( "or Z2 = " ); + pnum( &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" ); + Ptemp = Q / LOG(Radix); + pnum( &Ptemp ); + printf( "extra B-digits, i.e.\n" ); + Ptemp = Q / LOG(10.); + printf("roughly " ); + pnum( &Ptemp ); + printf( "extra significant decimals.\n"); + } + 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\n\ +or 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(" subtration 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)) { + printf("Y2 = "); + pnum( &Y2 ); + printf("Y1 = "); + pnum( &Y1 ); + printf("U2 = "); + pnum( &U2 ); + 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); + Ptemp = 3.0; + Random9 = SQRT(Ptemp); + 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 < 10) { + 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 " ); + pnum( &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; + 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) { |