summaryrefslogtreecommitdiff
path: root/libm/ldouble/lparanoi.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/ldouble/lparanoi.c')
-rw-r--r--libm/ldouble/lparanoi.c2348
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) {