diff options
Diffstat (limited to 'libm/double/paranoia.c')
| -rw-r--r-- | libm/double/paranoia.c | 2156 |
1 files changed, 0 insertions, 2156 deletions
diff --git a/libm/double/paranoia.c b/libm/double/paranoia.c deleted file mode 100644 index 49ff72623..000000000 --- a/libm/double/paranoia.c +++ /dev/null @@ -1,2156 +0,0 @@ -/* 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> - -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 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; -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; -volatile FLOAT VV; -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(); -} - -main() -{ - /* Set coprocessor to double precision, no arith traps. */ - /* __setfpucw(0x127f);*/ - dprec(); - /* 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 = %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."); - - /*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\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)) { - 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 < 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 %.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; - 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 |