forked from Imagelibrary/rtems
2310 lines
70 KiB
C
2310 lines
70 KiB
C
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#include <bsp.h>
|
|
/*
|
|
* 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).
|
|
*
|
|
* Sun May 16 18:21:51 MDT 1993 Jeffrey Wheat (cassidy@cygnus.com)
|
|
* Removed KR_headers defines, removed ANSI prototyping
|
|
* Cleaned up and reformated code. Added special CYGNUS
|
|
* "verbose" mode type messages (-DCYGNUS).
|
|
* Note: This code is VERY NASTY.
|
|
*
|
|
* Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
|
|
* compile with -DKR_headers or insert
|
|
* #define KR_headers
|
|
* at the beginning if you have an old-style C compiler.
|
|
*
|
|
* (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_PRECISION 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
|
|
*
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
|
|
/*
|
|
* To compile this on host using only libm from newlib (and using host libc)
|
|
* do:
|
|
* gcc -g -DNEED_REENT -DCYGNUS paranoia.c .../newlib-1.6/newlib/libm.a
|
|
*/
|
|
|
|
#ifdef NEED_REENT
|
|
#include <reent.h>
|
|
struct _reent libm_reent = _REENT_INIT(libm_reent);
|
|
struct _reent *_impure_ptr = &libm_reent;
|
|
#endif
|
|
|
|
#ifndef NOSIGNAL
|
|
#include <signal.h>
|
|
#include <setjmp.h>
|
|
#else /* NOSIGNAL */
|
|
#define longjmp(e,v)
|
|
#define setjmp(e) 0
|
|
#define jmp_buf int
|
|
#endif /* NOSIGNAL */
|
|
|
|
#ifdef SINGLE_PRECISION
|
|
#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 /* !SINGLE_PRECISION */
|
|
#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 /* SINGLE_PRECISION */
|
|
|
|
jmp_buf ovfl_buf;
|
|
/* extern double fabs (), floor (), log (), pow (), sqrt (); */
|
|
/* extern void exit (); */
|
|
extern void _sigfpe(int);
|
|
typedef void (*Sig_type) (int);
|
|
FLOAT Sign (FLOAT), Random (void);
|
|
extern void BadCond (int, char*);
|
|
extern void SqXMinX (int);
|
|
extern void TstCond (int, int, char*);
|
|
extern void notify (char *);
|
|
/* extern int read (); */
|
|
extern void Characteristics (void);
|
|
extern void Heading (void);
|
|
extern void History (void);
|
|
extern void Instructions (void);
|
|
extern void IsYeqX (void);
|
|
extern void NewD (void);
|
|
extern void Pause (void);
|
|
extern void PrintIfNPositive (void);
|
|
extern void SR3750 (void);
|
|
extern void SR3980 (void);
|
|
extern void TstPtUf (void);
|
|
|
|
void msglist(char**);
|
|
|
|
Sig_type sigsave;
|
|
|
|
#define KEYBOARD 0
|
|
|
|
FLOAT Radix, BInvrse, RadixD2, BMinusU2;
|
|
|
|
/*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 HVar, 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 WVar;
|
|
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
|
|
*/
|
|
|
|
int batchmode; /* global batchmode test */
|
|
|
|
/* program name and version variables and macro */
|
|
char *temp;
|
|
char *program_name;
|
|
char *program_vers;
|
|
|
|
#ifndef VERSION
|
|
#define VERSION "1.1 [cygnus]"
|
|
#endif /* VERSION */
|
|
|
|
#define basename(cp) ((temp=(char *)strrchr((cp), '/')) ? temp+1 : (cp))
|
|
|
|
#ifndef BATCHMODE
|
|
# ifdef CYGNUS
|
|
# define BATCHMODE
|
|
# endif
|
|
#endif
|
|
|
|
/* floating point exception receiver */
|
|
void
|
|
_sigfpe (x)
|
|
int x;
|
|
{
|
|
fpecount++;
|
|
printf ("\n* * * FLOATING-POINT ERROR %d * * *\n", x);
|
|
fflush (stdout);
|
|
if (sigsave) {
|
|
#ifndef NOSIGNAL
|
|
signal (SIGFPE, sigsave);
|
|
#endif /* NOSIGNAL */
|
|
sigsave = 0;
|
|
longjmp (ovfl_buf, 1);
|
|
}
|
|
exit (1);
|
|
}
|
|
|
|
#ifdef NOMAIN
|
|
#define main paranoia
|
|
int paranoia(int, char**);
|
|
#endif
|
|
|
|
int
|
|
main (argc, argv)
|
|
int argc;
|
|
char **argv;
|
|
{
|
|
/* 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;
|
|
|
|
#ifdef BATCHMODE
|
|
batchmode = 1; /* run test in batchmode? */
|
|
#else /* !BATCHMODE */
|
|
batchmode = 0; /* run test interactively */
|
|
#endif /* BATCHMODE */
|
|
|
|
program_name = basename (argv[0]);
|
|
program_vers = VERSION;
|
|
|
|
printf ("%s version %s\n", program_name, program_vers);
|
|
|
|
/*=============================================*/
|
|
Milestone = 0;
|
|
/*=============================================*/
|
|
#ifndef NOSIGNAL
|
|
signal (SIGFPE, _sigfpe);
|
|
#endif
|
|
|
|
if (!batchmode) {
|
|
Instructions ();
|
|
Pause ();
|
|
Heading ();
|
|
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");
|
|
/*=============================================*/
|
|
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");
|
|
WVar = One;
|
|
do {
|
|
WVar = WVar + WVar;
|
|
Y = WVar + One;
|
|
Z = Y - WVar;
|
|
Y = Z - One;
|
|
}
|
|
while (MinusOne + FABS (Y) < Zero);
|
|
/*.. now WVar is just big enough that |((WVar+1)-WVar)-1| >= 1 ...*/
|
|
Precision = Zero;
|
|
Y = One;
|
|
do {
|
|
Radix = WVar + Y;
|
|
Y = Y + Y;
|
|
Radix = Radix - WVar;
|
|
}
|
|
while (Radix == Zero);
|
|
if (Radix < Two)
|
|
Radix = One;
|
|
printf ("Radix = %f .\n", Radix);
|
|
if (Radix != 1) {
|
|
WVar = One;
|
|
do {
|
|
Precision = Precision + One;
|
|
WVar = WVar * Radix;
|
|
Y = WVar + One;
|
|
}
|
|
while ((Y - WVar) == One);
|
|
}
|
|
/*... now WVar == Radix^Precision is barely too big to satisfy (WVar+1)-WVar == 1
|
|
...*/
|
|
U1 = One / WVar;
|
|
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);
|
|
WVar = 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 ();
|
|
/*=============================================*/
|
|
Milestone = 35;
|
|
/*=============================================*/
|
|
if (Radix >= Two) {
|
|
X = WVar / (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 (" 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");
|
|
/*=============================================*/
|
|
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);
|
|
}
|
|
/*=============================================*/
|
|
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 != WVar - 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 != WVar - 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", WVar);
|
|
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;
|
|
Break = False;
|
|
do {
|
|
X = One;
|
|
SR3980 ();
|
|
if (I <= 10) {
|
|
I = 1023;
|
|
SR3980 ();
|
|
}
|
|
if (Z == MinusOne)
|
|
Break = True;
|
|
else {
|
|
Z = MinusOne;
|
|
/* .. if(-1)^N is invalid, replace MinusOne by One. */
|
|
I = -4;
|
|
}
|
|
}
|
|
while (!Break);
|
|
PrintIfNPositive ();
|
|
N1 = N;
|
|
N = 0;
|
|
Z = A1;
|
|
M = (int) FLOOR (Two * LOG (WVar) / 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 discrepancies found.\n");
|
|
if (N > 0)
|
|
Pause ();
|
|
else
|
|
printf ("\n");
|
|
/*=============================================*/
|
|
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;
|
|
HVar = One / HInvrse;
|
|
/* ... 1/HInvrse == HVar == Min(1/Radix, 1/2) */
|
|
CInvrse = One / C;
|
|
E0 = C;
|
|
Z = E0 * HVar;
|
|
/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
|
|
do {
|
|
Y = E0;
|
|
E0 = Z;
|
|
Z = E0 * HVar;
|
|
}
|
|
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 * HVar;
|
|
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 * HVar;
|
|
}
|
|
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 = HVar * 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 = HVar + HVar;
|
|
} else
|
|
R = SQRT (Underflow / UfThold);
|
|
sigsave = 0;
|
|
if (R <= HVar) {
|
|
Z = R * UfThold;
|
|
X = Z * (One + R * HVar * (One + HVar));
|
|
} else {
|
|
Z = UfThold;
|
|
X = Z * (One + HVar * HVar * (One + HVar));
|
|
}
|
|
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);
|
|
}
|
|
/*=============================================*/
|
|
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, Y2);
|
|
printf ("actually calculating yields:");
|
|
if (setjmp (ovfl_buf)) {
|
|
sigsave = 0;
|
|
BadCond (Serious, "trap on underflow.\n");
|
|
} else {
|
|
sigsave = _sigfpe;
|
|
V9 = POW (HInvrse, Y2);
|
|
sigsave = 0;
|
|
printf (" %.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 < WVar)
|
|
BadCond (Serious, "");
|
|
else
|
|
BadCond (Defect, "");
|
|
printf ("Comparison alleges that Z = %17e\n", Z);
|
|
printf (" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
|
|
}
|
|
}
|
|
/*=============================================*/
|
|
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 BATCHMODE
|
|
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 /* !BATCHMODE */
|
|
sigsave = _sigfpe;
|
|
printf (" Trying to compute 1 / 0 produces ...");
|
|
if (!setjmp (ovfl_buf))
|
|
printf (" %.7e .\n", One / MyZero);
|
|
sigsave = 0;
|
|
#ifndef BATCHMODE
|
|
} 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 /* !BATCHMODE */
|
|
sigsave = _sigfpe;
|
|
printf ("\n Trying to compute 0 / 0 produces ...");
|
|
if (!setjmp (ovfl_buf))
|
|
printf (" %.7e .\n", Zero / MyZero);
|
|
sigsave = 0;
|
|
#ifndef BATCHMODE
|
|
} else
|
|
printf ("O.K.\n");
|
|
#endif /* !BATCHMODE */
|
|
/*=============================================*/
|
|
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;
|
|
}
|
|
|
|
FLOAT
|
|
Sign (X)
|
|
FLOAT X;
|
|
{
|
|
return X >= 0. ? 1.0 : -1.0;
|
|
}
|
|
|
|
void
|
|
Pause ()
|
|
{
|
|
#ifndef BATCHMODE
|
|
char ch[8];
|
|
|
|
printf ("\nTo continue, press RETURN");
|
|
fflush (stdout);
|
|
read (KEYBOARD, ch, 8);
|
|
#endif /* !BATCHMODE */
|
|
#ifndef CYGNUS
|
|
printf ("\nDiagnosis resumes after milestone Number %d", Milestone);
|
|
printf (" Page: %d\n\n", PageNo);
|
|
++Milestone;
|
|
++PageNo;
|
|
#endif /* !CYGNUS */
|
|
}
|
|
|
|
void
|
|
TstCond (K, Valid, T)
|
|
int K, Valid;
|
|
char *T;
|
|
{
|
|
#ifdef CYGNUS
|
|
printf ("TEST: %s\n", T);
|
|
#endif /* CYGNUS */
|
|
if (!Valid) {
|
|
BadCond (K, T);
|
|
printf (".\n");
|
|
}
|
|
#ifdef CYGNUS
|
|
printf ("PASS: %s\n", T);
|
|
#endif /* CYGNUS */
|
|
}
|
|
|
|
void
|
|
BadCond (K, T)
|
|
int K;
|
|
char *T;
|
|
{
|
|
static char *msg[] =
|
|
{"FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW"};
|
|
|
|
ErrCnt[K] = ErrCnt[K] + 1;
|
|
#ifndef CYGNUS
|
|
printf ("%s: %s", msg[K], T);
|
|
#else
|
|
printf ("ERROR: Severity: %s: %s", msg[K], T);
|
|
#endif /* CYGNUS */
|
|
}
|
|
|
|
/*
|
|
* 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);
|
|
}
|
|
|
|
void
|
|
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");
|
|
}
|
|
}
|
|
|
|
void
|
|
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;
|
|
}
|
|
|
|
void
|
|
SR3750 ()
|
|
{
|
|
if (!((X - Radix < Z2 - Radix) || (X - Z2 > WVar - 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;
|
|
}
|
|
}
|
|
|
|
void
|
|
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. */
|
|
}
|
|
}
|
|
|
|
void
|
|
SR3980 ()
|
|
{
|
|
do {
|
|
Q = (FLOAT) I;
|
|
Y = POW (Z, Q);
|
|
IsYeqX ();
|
|
if (++I > M)
|
|
break;
|
|
X = Z * X;
|
|
}
|
|
while (X < WVar);
|
|
}
|
|
|
|
void
|
|
PrintIfNPositive ()
|
|
{
|
|
if (N > 0)
|
|
printf ("Similar discrepancies have occurred %d times.\n", N);
|
|
}
|
|
|
|
void
|
|
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 ();
|
|
}
|
|
}
|
|
}
|
|
|
|
void
|
|
notify (s)
|
|
char *s;
|
|
{
|
|
printf ("%s test appears to be inconsistent...\n", s);
|
|
printf (" PLEASE NOTIFY KARPINKSI!\n");
|
|
}
|
|
|
|
void
|
|
msglist (s)
|
|
char **s;
|
|
{
|
|
while (*s)
|
|
printf ("%s\n", *s++);
|
|
}
|
|
|
|
void
|
|
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);
|
|
}
|
|
|
|
void
|
|
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_PRECISION
|
|
"\tPrecision:\tsingle;",
|
|
#else /* !SINGLE_PRECISION */
|
|
"\tPrecision:\tdouble;",
|
|
#endif /* SINGLE_PRECISION */
|
|
"\tVersion:\t10 February 1989;",
|
|
"\tComputer:\n",
|
|
"\tCompiler:\n",
|
|
"\tOptimization level:\n",
|
|
"\tOther relevant compiler options:",
|
|
0};
|
|
|
|
msglist (head);
|
|
}
|
|
|
|
void
|
|
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);
|
|
}
|
|
|
|
void
|
|
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);
|
|
}
|