/*	A C version of Kahan's Floating Point Test "Paranoia"

			Thos Sumner, UCSF, Feb. 1985
			David Gay, BTL, Jan. 1986

	This is a rewrite from the Pascal version by

			B. A. Wichmann, 18 Jan. 1985

	(and does NOT exhibit good C programming style).

(C) Apr 19 1983 in BASIC version by:
	Professor W. M. Kahan,
	567 Evans Hall
	Electrical Engineering & Computer Science Dept.
	University of California
	Berkeley, California 94720
	USA

converted to Pascal by:
	B. A. Wichmann
	National Physical Laboratory
	Teddington Middx
	TW11 OLW
	UK

converted to C by:

	David M. Gay		and	Thos Sumner
	AT&T Bell Labs			Computer Center, Rm. U-76
	600 Mountain Avenue		University of California
	Murray Hill, NJ 07974		San Francisco, CA 94143
	USA				USA

with simultaneous corrections to the Pascal source (reflected
in the Pascal source available over netlib).

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/^GENERIC/extern &/
1,$s! = .*!;!
/^Guard/,/^Round/s/^/extern /
/^jmp_buf/s/^/extern /
/^Sig_type/s/^/extern /
a
extern void Sigfpe();
.
w paranoia.h
q

*/

#include "ucbtest.h"

#define NOSIGNAL
#define NOPAUSE

#include <stdio.h>
#ifndef NOSIGNAL
#include <signal.h>
#endif
#include <setjmp.h>
#if defined(__STDC__) || (sunos==4) || (sunos==5)
#include <stdlib.h>
#endif

jmp_buf ovfl_buf;
typedef void (*Sig_type)();
Sig_type sigsave;

#define KEYBOARD 0

GENERIC Radix, BInvrse, RadixD2, BMinusU2;
GENERIC Sign(), Random();

/*Small floating point constants.*/
GENERIC Zero = 0.0;
GENERIC Half = 0.5;
GENERIC One = 1.0;
GENERIC Two = 2.0;
GENERIC Three = 3.0;
GENERIC Four = 4.0;
GENERIC Five = 5.0;
GENERIC Eight = 8.0;
GENERIC Nine = 9.0;
GENERIC TwentySeven = 27.0;
GENERIC ThirtyTwo = 32.0;
GENERIC TwoForty = 240.0;
GENERIC MinusOne = -1.0;
GENERIC 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];
GENERIC AInvrse, A1;
GENERIC C, CInvrse;
GENERIC D, FourD;
GENERIC E0, E1, Exp2, E3, MinSqEr;
GENERIC SqEr, MaxSqEr, E9;
GENERIC Third;
GENERIC F6, F9;
GENERIC H, HInvrse;
int I;
GENERIC StickyBit, J;
GENERIC MyZero;
GENERIC Precision;
GENERIC Q, Q9;
GENERIC R, Random9;
GENERIC T, Underflow, S;
GENERIC OneUlp, UfThold, U1, U2;
GENERIC V, V0, V9;
GENERIC W;
GENERIC X, X1, X2, X8, Random1;
GENERIC Y, Y1, Y2, Random2;
GENERIC Z, PseudoZero, Z1, Z2, Z9;
int ErrCnt[4];
int fpecount;
int Milestone;
int PageNo;
int M, N, N1;
Guard GMult, GDiv, GAddSub;
Rounding RMult, RDiv, RAddSub, RSqrt;
int Break, Done, NotMonot, Monot, Anomaly, IEEE,
		SqRWrng, UfNGrad;
/* Computed constants. */
/*U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
/*U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */

/* floating point exception receiver */
void
Sigfpe(arg)
	int arg;
{
	fpecount++;
	printf("\n* * * GENERICING-POINT ERROR * * *\n");
	fflush(stdout);
	if (sigsave) {
#ifndef NOSIGNAL
		signal(SIGFPE, sigsave);
#endif
		sigsave = 0;
		longjmp(ovfl_buf, 1);
		}
	abort();
}

void Instructions(), Pause(), Heading(), Characteristics(), History();
void TstCond(), TstPtUf(), BadCond(), notify(), SqXMinX(), NewD();
void SR3750(), SR3980(), PrintIfNPositive(), IsYeqX(), msglist();

#ifdef __STDC__
extern int read(int, char *, unsigned);
#endif /* __STDC__ */

main()
{
	(void) ucbstart( __FILE__ , __LINE__ );

	/* First two assignments use integer right-hand sides. */
	
	Zero = 0;
	One = 1;
	Two = One + One;
	Three = Two + One;
	Four = Three + One;
	Five = Four + One;
	Eight = Four + Four;
	Nine = Three * Three;
	TwentySeven = Nine * Three;
	ThirtyTwo = Four * Eight;
	TwoForty = Four * Five * Three * Four;
	MinusOne = -One;
	Half = One / Two;
	OneAndHalf = One + Half;
	ErrCnt[Failure] = 0;
	ErrCnt[Serious] = 0;
	ErrCnt[Defect] = 0;
	ErrCnt[Flaw] = 0;
	PageNo = 1;
	/*=============================================*/
	Milestone = 0;
	/*=============================================*/
#ifndef NOSIGNAL
	signal(SIGFPE, Sigfpe);
#endif
	Instructions();
	Pause();
	Heading();
	Pause();
	Characteristics();
	Pause();
	History();
	Pause();
	/*=============================================*/
	Milestone = 7;
	/*=============================================*/
	printf("Program is now RUNNING tests on small integers:\n");
	
	TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
		   && (One > Zero) && (One + One == Two),
			"0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
	Z = - Zero;
	if (Z != 0.0) {
		ErrCnt[Failure] = ErrCnt[Failure] + 1;
		printf("Comparison alleges that -0.0 is Non-zero!\n");
		U2 = 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 + ABS(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();
#if (sunpro_version >= 100) || ((sunos_version/100) == 4) 
	ieee_retrospective_();
#endif	
	}
#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 + ABS(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;
#ifdef QP
	printf("Radix = %Lf .\n", Radix);
#else
	printf("Radix = %f .\n", Radix);
#endif
	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;
#ifdef QP
	printf("Closest relative separation found is U1 = %.7Le .\n\n", U1);
#else
	printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
#endif
	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 = ABS(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 = ABS(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");
# ifdef QP
	else printf("gets better closest relative separation U1 = %.7Le .\n", U1);
# else
	else printf("gets better closest relative separation U1 = %.7e .\n", U1);
# endif
	W = One / U1;
	F9 = (Half - U1) + Half;
	Radix = FLOOR(0.01 + U2 / U1);
	if (Radix == E0) printf("Radix confirmed.\n");
# ifdef QP
	else printf("MYSTERY: recalculated Radix = %.7Le .\n", Radix);
# else
	else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
# endif
	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 (ABS(X - Y) * Four < One) X = Y;
		Precision = X / TwoForty;
		Y = FLOOR(Half + Precision);
		if (ABS(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");
#ifdef QP
	else printf("The number of significant digits of the Radix is %Lf .\n",
			Precision);
#else
	else printf("The number of significant digits of the Radix is %f .\n",
			Precision);
#endif
	TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
		   "Precision worse than 5 decimal figures  ");
	/*=============================================*/
	Milestone = 30;
	/*=============================================*/
	/* Test for extra-precise subepressions */
	X = ABS(((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 = ABS((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");
# ifdef QP
		printf("respectively  %.7Le,  %.7Le,  %.7Le,\n", X1, Y1, Z1);
# else
		printf("respectively  %.7e,  %.7e,  %.7e,\n", X1, Y1, Z1);
# endif
		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");
# ifdef QP
				printf("\tU1 = %.7Le, Z1 - U1 = %.7Le\n",U1,Z1-U1);
				printf("\tU2 = %.7Le, Z2 - U2 = %.7Le\n",U2,Z2-U2);
# else
				printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
				printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
# endif
				}
			else {
				if ((Z1 <= Zero) || (Z2 <= Zero)) {
# ifdef QP
					printf("Because of unusual Radix = %Lf", Radix);
					printf(", or exact rational arithmetic a result\n");
					printf("Z1 = %.7Le, or Z2 = %.7Le ", Z1, Z2);
# else
					printf("Because of unusual Radix = %f", Radix);
					printf(", or exact rational arithmetic a result\n");
					printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
# endif
					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");
# ifdef QP
					printf("precisely with about %Lg extra B-digits, i.e.\n",
						(Q / LOG(Radix)));
					printf("roughly %Lg extra significant decimals.\n",
						Q / LOG(10.));
# else
					printf("precisely with about %g extra B-digits, i.e.\n",
						(Q / LOG(Radix)));
					printf("roughly %g extra significant decimals.\n",
						Q / LOG(10.));
# endif
					}
				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 = ABS((X + Z) - X * X) - U2;
	X = One - U2;
	Z = ABS((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");
/* dgh */	RSqrt = Other ;
	RMult = Other;
	RDiv = Other;
	RAddSub = Other;
	RadixD2 = Radix / Two;
	A1 = Two;
	Done = False;
	do  {
		AInvrse = Radix;
		do  {
			X = AInvrse;
			AInvrse = AInvrse / A1;
			} while ( ! (FLOOR(AInvrse) != AInvrse));
		Done = (X == One) || (A1 > Three);
		if (! Done) A1 = Nine + One;
		} while ( ! (Done));
	if (X == One) A1 = Radix;
	AInvrse = One / A1;
	X = A1;
	Y = AInvrse;
	Done = False;
	do  {
		Z = X * Y - Half;
		TstCond (Failure, Z == Half,
			"X * (1/X) differs from 1");
		Done = X == Radix;
		X = Radix;
		Y = One / X;
		} while ( ! (Done));
	Y2 = One + U2;
	Y1 = One - U2;
	X = OneAndHalf - U2;
	Y = OneAndHalf + U2;
	Z = (X - U2) * Y2;
	T = Y * Y1;
	Z = Z - X;
	T = T - X;
	X = X * Y2;
	Y = (Y + U2) * Y1;
	X = X - OneAndHalf;
	Y = Y - OneAndHalf;
	if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
		X = (OneAndHalf + U2) * Y2;
		Y = OneAndHalf - U2 - U2;
		Z = OneAndHalf + U2 + U2;
		T = (OneAndHalf - U2) * Y1;
		X = X - (Z + U2);
		StickyBit = Y * Y1;
		S = Z * Y2;
		T = T - Y;
		Y = (U2 - Y) + StickyBit;
		Z = S - (Z + U2 + U2);
		StickyBit = (Y2 + U2) * Y1;
		Y1 = Y2 * Y1;
		StickyBit = StickyBit - Y2;
		Y1 = Y1 - Half;
		if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
			&& ( StickyBit == Zero) && (Y1 == Half)) {
			RMult = Rounded;
			printf("Multiplication appears to round correctly.\n");
			}
		else	if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
				&& (T < Zero) && (StickyBit + U2 == Zero)
				&& (Y1 < Half)) {
				RMult = Chopped;
				printf("Multiplication appears to chop.\n");
				}
			else printf("* is neither chopped nor correctly rounded.\n");
		if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
		}
	else printf("* is neither chopped nor correctly rounded.\n");
	/*=============================================*/
	Milestone = 45;
	/*=============================================*/
	Y2 = One + U2;
	Y1 = One - U2;
	Z = OneAndHalf + U2 + U2;
	X = Z / Y2;
	T = OneAndHalf - U2 - U2;
	Y = (T - U2) / Y1;
	Z = (Z + U2) / Y2;
	X = X - OneAndHalf;
	Y = Y - T;
	T = T / Y1;
	Z = Z - (OneAndHalf + U2);
	T = (U2 - OneAndHalf) + T;
	if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
		X = OneAndHalf / Y2;
		Y = OneAndHalf - U2;
		Z = OneAndHalf + U2;
		X = X - Y;
		T = OneAndHalf / Y1;
		Y = Y / Y1;
		T = T - (Z + U2);
		Y = Y - Z;
		Z = Z / Y2;
		Y1 = (Y2 + U2) / Y2;
		Z = Z - OneAndHalf;
		Y2 = Y1 - Y2;
		Y1 = (F9 - U1) / F9;
		if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
			&& (Y2 == Zero) && (Y2 == Zero)
			&& (Y1 - Half == F9 - Half )) {
			RDiv = Rounded;
			printf("Division appears to round correctly.\n");
			if (GDiv == No) notify("Division");
			}
		else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
			&& (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
			RDiv = Chopped;
			printf("Division appears to chop.\n");
			}
		}
	if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
	BInvrse = One / Radix;
	TstCond (Failure, (BInvrse * Radix - Half == Half),
		   "Radix * ( 1 / Radix ) differs from 1");
	/*=============================================*/
	/*SPLIT
	}
#include "paranoia.h"
part4(){
*/
	Milestone = 50;
	/*=============================================*/
	TstCond (Failure, ((F9 + U1) - Half == Half)
		   && ((BMinusU2 + U2 ) - One == Radix - One),
		   "Incomplete carry-propagation in Addition");
	X = One - U1 * U1;
	Y = One + U2 * (One - U2);
	Z = F9 - Half;
	X = (X - Half) - Z;
	Y = Y - One;
	if ((X == Zero) && (Y == Zero)) {
		RAddSub = Chopped;
		printf("Add/Subtract appears to be chopped.\n");
		}
	if (GAddSub == Yes) {
		X = (Half + U2) * U2;
		Y = (Half - U2) * U2;
		X = One + X;
		Y = One + Y;
		X = (One + U2) - X;
		Y = One - Y;
		if ((X == Zero) && (Y == Zero)) {
			X = (Half + U2) * U1;
			Y = (Half - U2) * U1;
			X = One - X;
			Y = One - Y;
			X = F9 - X;
			Y = One - Y;
			if ((X == Zero) && (Y == Zero)) {
				RAddSub = Rounded;
				printf("Addition/Subtraction appears to round correctly.\n");
				if (GAddSub == No) notify("Add/Subtract");
				}
			else printf("Addition/Subtraction neither rounds nor chops.\n");
			}
		else printf("Addition/Subtraction neither rounds nor chops.\n");
		}
	else printf("Addition/Subtraction neither rounds nor chops.\n");
	S = One;
	X = One + Half * (One + Half);
	Y = (One + U2) * Half;
	Z = X - Y;
	T = Y - X;
	StickyBit = Z + T;
	if (StickyBit != Zero) {
		S = Zero;
		BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
		}
	StickyBit = Zero;
	if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
		&& (RMult == Rounded) && (RDiv == Rounded)
		&& (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
		printf("Checking for sticky bit.\n");
		X = (Half + U1) * U2;
		Y = Half * U2;
		Z = One + Y;
		T = One + X;
		if ((Z - One <= Zero) && (T - One >= U2)) {
			Z = T + Y;
			Y = Z - X;
			if ((Z - T >= U2) && (Y - T == Zero)) {
				X = (Half + U1) * U1;
				Y = Half * U1;
				Z = One - Y;
				T = One - X;
				if ((Z - One == Zero) && (T - F9 == Zero)) {
					Z = (Half - U1) * U1;
					T = F9 - Z;
					Q = F9 - Y;
					if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
						Z = (One + U2) * OneAndHalf;
						T = (OneAndHalf + U2) - Z + U2;
						X = One + Half / Radix;
						Y = One + Radix * U2;
						Z = X * Y;
						if (T == Zero && X + Radix * U2 - Z == Zero) {
							if (Radix != Two) {
								X = Two + U2;
								Y = X / Two;
								if ((Y - One == Zero)) StickyBit = S;
								}
							else StickyBit = S;
							}
						}
					}
				}
			}
		}
	if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
	else printf("Sticky bit used incorrectly or not at all.\n");
	TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
			RMult == Other || RDiv == Other || RAddSub == Other),
		"lack(s) of guard digits or failure(s) to correctly round or chop\n\
(noted above) count as one flaw in the final tally below");
	/*=============================================*/
	Milestone = 60;
	/*=============================================*/
	printf("\n");
	printf("Does Multiplication commute?  ");
	printf("Testing on %d random pairs.\n", NoTrials);
	Random9 = SQRT(3.0);
	Random1 = Third;
	I = 1;
	do  {
		X = Random();
		Y = Random();
		Z9 = Y * X;
		Z = X * Y;
		Z9 = Z - Z9;
		I = I + 1;
		} while ( ! ((I > NoTrials) || (Z9 != Zero)));
	if (I == NoTrials) {
		Random1 = One + Half / Three;
		Random2 = (U2 + U1) + One;
		Z = Random1 * Random2;
		Y = Random2 * Random1;
		Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
			Three) * ((U2 + U1) + One);
		}
	if (! ((I == NoTrials) || (Z9 == Zero)))
		BadCond(Defect, "X * Y == Y * X trial fails.\n");
	else printf("     No failures found in %d integer pairs.\n", NoTrials);
	/*=============================================*/
	Milestone = 70;
	/*=============================================*/
	printf("\nRunning test of square root(x).\n");
	TstCond (Failure, (Zero == SQRT(Zero))
		   && (- Zero == SQRT(- Zero))
		   && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
	MinSqEr = Zero;
	MaxSqEr = Zero;
	J = Zero;
	X = Radix;
	OneUlp = U2;
	SqXMinX (Serious);
	X = BInvrse;
	OneUlp = BInvrse * U1;
	SqXMinX (Serious);
	X = U1;
	OneUlp = U1 * U1;
	SqXMinX (Serious);
	if (J != Zero) Pause();
	printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
	J = Zero;
	X = Two;
	Y = Radix;
	if ((Radix != One)) do  {
		X = Y;
		Y = Radix * Y;
		} while ( ! ((Y - X >= NoTrials)));
	OneUlp = X * U2;
	I = 1;
	while (I < 10) {
		X = X + One;
		SqXMinX (Defect);
		if (J > Zero) break;
		I = I + 1;
		}
	printf("Test for sqrt monotonicity.\n");
	I = - 1;
	X = BMinusU2;
	Y = Radix;
	Z = Radix + Radix * U2;
	NotMonot = False;
	Monot = False;
	while ( ! (NotMonot || Monot)) {
		I = I + 1;
		X = SQRT(X);
		Q = SQRT(Y);
		Z = SQRT(Z);
		if ((X > Q) || (Q > Z)) NotMonot = True;
		else {
			Q = FLOOR(Q + Half);
			if ((I > 0) || (Radix == Q * Q)) Monot = True;
			else if (I > 0) {
			if (I > 1) Monot = True;
			else {
				Y = Y * BInvrse;
				X = Y - U1;
				Z = Y + U1;
				}
			}
			else {
				Y = Q;
				X = Y - U2;
				Z = Y + U2;
				}
			}
		}
	if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
	else {
		BadCond(Defect, "");
# ifdef QP
		printf("sqrt(X) is non-monotonic for X near %.7Le .\n", Y);
# else
		printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
# endif
		}
	/*=============================================*/
	/*SPLIT
	}
#include "paranoia.h"
part5(){
*/
	Milestone = 80;
	/*=============================================*/
	MinSqEr = MinSqEr + Half;
	MaxSqEr = MaxSqEr - Half;
	Y = (SQRT(One + U2) - One) / U2;
	SqEr = (Y - One) + U2 / Eight;
	if (SqEr > MaxSqEr) MaxSqEr = SqEr;
	SqEr = Y + U2 / Eight;
	if (SqEr < MinSqEr) MinSqEr = SqEr;
	Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
	SqEr = Y + U1 / Eight;
	if (SqEr > MaxSqEr) MaxSqEr = SqEr;
	SqEr = (Y + One) + U1 / Eight;
	if (SqEr < MinSqEr) MinSqEr = SqEr;
	OneUlp = U2;
	X = OneUlp;
	for( Indx = 1; Indx <= 3; ++Indx) {
		Y = SQRT((X + U1 + X) + F9);
		Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
		Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
		SqEr = (Y + Half) + Z;
		if (SqEr < MinSqEr) MinSqEr = SqEr;
		SqEr = (Y - Half) + Z;
		if (SqEr > MaxSqEr) MaxSqEr = SqEr;
		if (((Indx == 1) || (Indx == 3))) 
			X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
		else {
			OneUlp = U1;
			X = - OneUlp;
			}
		}
	/*=============================================*/
	Milestone = 85;
	/*=============================================*/
	SqRWrng = False;
	Anomaly = False;
	if (Radix != One) {
		printf("Testing whether sqrt is rounded or chopped.\n");
		D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
	/* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
		X = D / Radix;
		Y = D / A1;
		if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
			Anomaly = True;
			}
		else {
			X = Zero;
			Z2 = X;
			Y = One;
			Y2 = Y;
			Z1 = Radix - One;
			FourD = Four * D;
			do  {
				if (Y2 > Z2) {
					Q = Radix;
					Y1 = Y;
					do  {
						X1 = ABS(Q + FLOOR(Half - Q / Y1) * Y1);
						Q = Y1;
						Y1 = X1;
						} while ( ! (X1 <= Zero));
					if (Q <= One) {
						Z2 = Y2;
						Z = Y;
						}
					}
				Y = Y + Two;
				X = X + Eight;
				Y2 = Y2 + X;
				if (Y2 >= FourD) Y2 = Y2 - FourD;
				} while ( ! (Y >= D));
			X8 = FourD - Z2;
			Q = (X8 + Z * Z) / FourD;
			X8 = X8 / Eight;
			if (Q != FLOOR(Q)) Anomaly = True;
			else {
				Break = False;
				do  {
					X = Z1 * Z;
					X = X - FLOOR(X / Radix) * Radix;
					if (X == One) 
						Break = True;
					else
						Z1 = Z1 - One;
					} while ( ! (Break || (Z1 <= Zero)));
				if ((Z1 <= Zero) && (! Break)) Anomaly = True;
				else {
					if (Z1 > RadixD2) Z1 = Z1 - Radix;
					do  {
						NewD();
						} while ( ! (U2 * D >= F9));
					if (D * Radix - D != W - D) Anomaly = True;
					else {
						Z2 = D;
						I = 0;
						Y = D + (One + Z) * Half;
						X = D + Z + Q;
						SR3750();
						Y = D + (One - Z) * Half + D;
						X = D - Z + D;
						X = X + Q + X;
						SR3750();
						NewD();
						if (D - Z2 != W - Z2) Anomaly = True;
						else {
							Y = (D - Z2) + (Z2 + (One - Z) * Half);
							X = (D - Z2) + (Z2 - Z + Q);
							SR3750();
							Y = (One + Z) * Half;
							X = Q;
							SR3750();
							if (I == 0) Anomaly = True;
							}
						}
					}
				}
			}
		if ((I == 0) || Anomaly) {
			BadCond(Failure, "Anomalous arithmetic with Integer < ");
# ifdef QP
			printf("Radix^Precision = %.7Le\n", W);
# else
			printf("Radix^Precision = %.7e\n", W);
# endif
			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");
# ifdef QP
		printf("Observed errors run from %.7Le ", MinSqEr - Half);
		printf("to %.7Le ulps.\n", Half + MaxSqEr);
# else
		printf("Observed errors run from %.7e ", MinSqEr - Half);
		printf("to %.7e ulps.\n", Half + MaxSqEr);
# endif
		TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
			"sqrt gets too many last digits wrong");
		}
	/*=============================================*/
	Milestone = 90;
	/*=============================================*/
	Pause();
	printf("Testing powers Z^i for small Integers Z and i.\n");
	N = 0;
	/* ... test powers of zero. */
	I = 0;
	Z = -Zero;
	M = 3.0;
	Break = False;
	do  {
		X = One;
		SR3980();
		if (I <= 10) {
			I = 1023;
			SR3980();
			}
		if (Z == MinusOne) Break = True;
		else {
			Z = MinusOne;
			PrintIfNPositive();
			N = 0;
			/* .. if(-1)^N is invalid, replace MinusOne by One. */
			I = - 4;
			}
		} while ( ! Break);
	PrintIfNPositive();
	N1 = N;
	N = 0;
	Z = A1;
	M = FLOOR(Two * LOG(W) / LOG(A1));
	Break = False;
	do  {
		X = Z;
		I = 1;
		SR3980();
		if (Z == AInvrse) Break = True;
		else Z = AInvrse;
		} while ( ! (Break));
	/*=============================================*/
		Milestone = 100;
	/*=============================================*/
	/*  Powers of Radix have been tested, */
	/*         next try a few primes     */
	M = NoTrials;
	Z = Three;
	do  {
		X = Z;
		I = 1;
		SR3980();
		do  {
			Z = Z + Two;
			} while ( Three * FLOOR(Z / Three) == Z );
		} while ( Z < Eight * Three );
	if (N > 0) {
		printf("Errors like this may invalidate financial calculations\n");
		printf("\tinvolving interest rates.\n");
		}
	PrintIfNPositive();
	N += N1;
	if (N == 0) printf("... no discrepancies found.\n");
	if (N > 0) Pause();
	else printf("\n");
	/*=============================================*/
	/*SPLIT
	}
#include "paranoia.h"
part6(){
*/
	Milestone = 110;
	/*=============================================*/
	printf("Seeking Underflow thresholds UfThold and E0.\n");
	D = U1;
	if (Precision != FLOOR(Precision)) {
		D = BInvrse;
		X = Precision;
		do  {
			D = D * BInvrse;
			X = X - One;
			} while ( X > Zero);
		}
	Y = One;
	Z = D;
	/* ... D is power of 1/Radix < 1. */
	do  {
		C = Y;
		Y = Z;
		Z = Y * Y;
		} while ((Y > Z) && (Z + Z > Z));
	Y = C;
	Z = Y * D;
	do  {
		C = Y;
		Y = Z;
		Z = Y * D;
		} while ((Y > Z) && (Z + Z > Z));
	if (Radix < Two) HInvrse = Two;
	else HInvrse = Radix;
	H = One / HInvrse;
	/* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
	CInvrse = One / C;
	E0 = C;
	Z = E0 * H;
	/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
	do  {
		Y = E0;
		E0 = Z;
		Z = E0 * H;
		} while ((E0 > Z) && (Z + Z > Z));
	UfThold = E0;
	E1 = Zero;
	Q = Zero;
	E9 = U2;
	S = One + E9;
	D = C * S;
	if (D <= C) {
		E9 = Radix * U2;
		S = One + E9;
		D = C * S;
		if (D <= C) {
			BadCond(Failure, "multiplication gets too many last digits wrong.\n");
			Underflow = E0;
			Y1 = Zero;
			PseudoZero = Z;
			Pause();
			}
		}
	else {
		Underflow = D;
		PseudoZero = Underflow * H;
		UfThold = Zero;
		do  {
			Y1 = Underflow;
			Underflow = PseudoZero;
			if (E1 + E1 <= E1) {
				Y2 = Underflow * HInvrse;
				E1 = ABS(Y1 - Y2);
				Q = Y1;
				if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
				}
			PseudoZero = PseudoZero * H;
			} while ((Underflow > PseudoZero)
				&& (PseudoZero + PseudoZero > PseudoZero));
		}
	/* Comment line 4530 .. 4560 */
	if (PseudoZero != Zero) {
		printf("\n");
		Z = PseudoZero;
	/* ... Test PseudoZero for "phoney- zero" violates */
	/* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
		   ... */
		if (PseudoZero <= Zero) {
			BadCond(Failure, "Positive expressions can underflow to an\n");
			printf("allegedly negative value\n");
# ifdef QP
			printf("PseudoZero that prints out as: %Lg .\n", PseudoZero);
# else
			printf("PseudoZero that prints out as: %g .\n", PseudoZero);
# endif
			X = - PseudoZero;
			if (X <= Zero) {
				printf("But -PseudoZero, which should be\n");
# ifdef QP
				printf("positive, isn't; it prints out as  %Lg .\n", X);
# else
				printf("positive, isn't; it prints out as  %g .\n", X);
# endif
				}
			}
		else {
			BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
# ifdef QP
			printf("value PseudoZero that prints out as %Lg .\n", PseudoZero);
# else
			printf("value PseudoZero that prints out as %g .\n", PseudoZero);
# endif
			}
		TstPtUf();
		}
	/*=============================================*/
	Milestone = 120;
	/*=============================================*/
	if (CInvrse * Y > CInvrse * Y1) {
		S = H * S;
		E0 = Underflow;
		}
	if (! ((E1 == Zero) || (E1 == E0))) {
		BadCond(Defect, "");
		if (E1 < E0) {
			printf("Products underflow at a higher");
			printf(" threshold than differences.\n");
			if (PseudoZero == Zero) 
			E0 = E1;
			}
		else {
			printf("Difference underflows at a higher");
			printf(" threshold than products.\n");
			}
		}
#ifdef QP
	printf("Smallest strictly positive number found is E0 = %Lg .\n", E0);
#else
	printf("Smallest strictly positive number found is E0 = %g .\n", E0);
#endif
	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");
# ifdef QP
			printf("approach a threshold = %.17Le\n", UfThold);;
			printf(" coming down from %.17Le\n", C);
# else
			printf("approach a threshold = %.17e\n", UfThold);;
			printf(" coming down from %.17e\n", C);
# endif
			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");
# ifdef QP
		printf("print out as Q = %.17Le, Y = %.17Le .\n", Q, Y);
		printf ("|Q - Y| = %.17Le .\n" , ABS(Q - Y2));
# else
		printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y);
		printf ("|Q - Y| = %.17e .\n" , ABS(Q - Y2));
# endif
		UfThold = Q;
		break;
	
		case	3:
		X = X;
		break;
	
		case	4:
		if ((Q == UfThold) && (E1 == E0)
			&& (ABS( UfThold - E1 / E9) <= E1)) {
#ifdef __STDC__
volatile
#endif
			GENERIC y1, y2, y3, y4;

			UfNGrad = False;
			printf("Underflow is gradual; it incurs Absolute Error =\n");
			printf("(roundoff in UfThold) < E0.\n");
			y1 = E0 * CInvrse;
			y2 = y1 * (OneAndHalf + U2);
			y3 = CInvrse * (One + U2);
			y4 = y2 / y3;
#ifdef QP
			printf(" test double rounding in gradual underflow E0 %Lg y4-E0 %Lg y %Lg %Lg %Lg %Lg \n",E0,y4-E0,y1,y2,y3,y4);
#else
			printf(" test double rounding in gradual underflow E0 %g y4-E0 %g y %g %g %g %g \n",E0,y4-E0,y1,y2,y3,y4);
#endif
			IEEE = (y4 == E0);
			}
		}
	if (UfNGrad) {
		printf("\n");
		R = SQRT(Underflow / UfThold);
		if (R <= H) {
			Z = R * UfThold;
			X = Z * (One + R * H * (One + H));
			}
		else {
			Z = UfThold;
			X = Z * (One + H * H * (One + H));
			}
		if (! ((X == Z) || (X - Z != Zero))) {
			BadCond(Flaw, "");
# ifdef QP
			printf("X = %.17Le\n\tis not equal to Z = %.17Le .\n", X, Z);
			Z9 = X - Z;
			printf("yet X - Z yields %.17Le .\n", Z9);
# else
			printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
			Z9 = X - Z;
			printf("yet X - Z yields %.17e .\n", Z9);
# endif
			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");
# ifdef QP
			printf("X / Z = 1 + %Lg .\n", (X / Z - Half) - Half);
# else
			printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
# endif
			}
		}
#ifdef QP
	printf("The Underflow threshold is %.17Le, %s\n", UfThold,
		   " below which");
#else
	printf("The Underflow threshold is %.17e, %s\n", UfThold,
		   " below which");
#endif
	printf("calculation may suffer larger Relative error than ");
	printf("merely roundoff.\n");
	Y2 = U1 * U1;
	Y = Y2 * Y2;
	Y2 = Y * U1;
	if (Y2 <= UfThold) {
		if (Y > E0) {
			BadCond(Defect, "");
			I = 5;
			}
		else {
			BadCond(Serious, "");
			I = 4;
			}
		printf("Range is too narrow; U1^%d Underflows.\n", I);
		}
	/*=============================================*/
	/*SPLIT
	}
#include "paranoia.h"
part7(){
*/
	Milestone = 130;
	/*=============================================*/
	Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
	Y2 = Y + Y;
	printf("Since underflow occurs below the threshold\n");
#ifdef QP
	printf("UfThold = (%.17Le) ^ (%.17Le)\nonly underflow ", HInvrse, Y);
	printf("should afflict the expression\n\t(%.17Le) ^ (%.17Le);\n", HInvrse, Y);
#else
	printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
	printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
#endif
	V9 = POW(HInvrse, Y2);

#ifdef QP
	printf("actually calculating yields: %.17Le .\n", V9);
#else
	printf("actually calculating yields: %.17e .\n", V9);
#endif
	if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
		BadCond(Serious, "this is not between 0 and underflow\n");
# ifdef QP
		printf("   threshold = %.17Le .\n", UfThold);
# else
		printf("   threshold = %.17e .\n", UfThold);
# endif
		}
	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");
# ifdef QP
		printf("   threshold = %.17Le .\n", UfThold);
# else
		printf("   threshold = %.17e .\n", UfThold);
# endif
		}
	/*=============================================*/
	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;
#ifdef QP
	printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17Le as X -> 1.\n",
		Exp2);
#else
	printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
		Exp2);
#endif
	for(I = 1;;) {
		Z = X - BInvrse;
		Z = (X + One) / (Z - (One - BInvrse));
		Q = POW(X, Z) - Exp2;
		if (ABS(Q) > TwoForty * U2) {
			N = 1;
	 		V9 = (X - BInvrse) - (One - BInvrse);
			BadCond(Defect, "Calculated");
# ifdef QP
			printf(" %.17Le for\n", POW(X,Z));
			printf("\t(1 + (%.17Le) ^ (%.17Le);\n", V9, Z);
			printf("\tdiffers from correct value by %.17Le .\n", Q);
# else
			printf(" %.17e for\n", POW(X,Z));
			printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
			printf("\tdiffers from correct value by %.17e .\n", Q);
# endif
			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");
	sigsave = Sigfpe;
	I = 0;
	Y = - CInvrse;
	V9 = HInvrse * Y;
	if (setjmp(ovfl_buf)) goto overflow;
	do {
		V = Y;
		Y = V9;
		V9 = HInvrse * Y;
		} while(V9 < Y);
	I = 1;
overflow:
	Z = V9;
	printf("Can `Z = -Y' overflow?\n");
#ifdef QP
	printf("Trying it on Y = %.17Le .\n", Y);
#else
	printf("Trying it on Y = %.17e .\n", Y);
#endif
	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, "");
# ifdef QP
		printf("overflow past %.17Le\n\tshrinks to %.17Le .\n", Y, Z);
# else
		printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
# endif
		}
	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;
#ifdef QP
	printf("Overflow threshold is V  = %.17Le .\n", V);
#else
	printf("Overflow threshold is V  = %.17e .\n", V);
#endif

	if (I) 
#ifdef QP
	        printf("Overflow saturates at V0 = %.17Le .\n", V0);
#else
	        printf("Overflow saturates at V0 = %.17e .\n", V0);
#endif
	else printf("There is no saturation value because \
the system traps on overflow.\n");

	V9 = V * One;

#ifdef QP
	printf("No Overflow should be signaled for V * 1 = %.17Le\n", V9);
#else
	printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
#endif

	V9 = V / One;

#ifdef QP
	printf("                           nor for V / 1 = %.17Le .\n", V9);
#else
	printf("                           nor for V / 1 = %.17e .\n", V9);
#endif

	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 ");
# ifdef QP
		printf("+-%Lg, +-%Lg\nand +-%Lg are confused by Overflow.",
			V, V0, UfThold);
# else
		printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
			V, V0, UfThold);
# endif
		}
	/*=============================================*/
	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) {
				if (V9 > U1) BadCond(Serious, "");
				else BadCond(Defect, "");
# ifdef QP
				printf("Comparison alleges that what prints as Z = %.17Le\n", Z);
				printf(" is too far from sqrt(Z) ^ 2 = %.17Le .\n", Y);
# else
				printf("Comparison alleges that what prints as Z = %.17e\n", Z);
				printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
# endif
				}
			}
		}
	/*=============================================*/
	Milestone = 180;
	/*=============================================*/
	for(Indx = 1; Indx <= 2; ++Indx) {
		if (Indx == 1) Z = V;
		else Z = V0;
		V9 = SQRT(Z);
		X = (One - Radix * E9) * V9;
		V9 = V9 * X;
		if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
			Y = V9;
			if (X < W) BadCond(Serious, "");
			else BadCond(Defect, "");
# ifdef QP
			printf("Comparison alleges that Z = %17Le\n", Z);
			printf(" is too far from sqrt(Z) ^ 2 (%.17Le) .\n", Y);
# else
			printf("Comparison alleges that Z = %17e\n", Z);
			printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
# endif
			}
		}
	/*=============================================*/
	/*SPLIT
	}
#include "paranoia.h"
part8(){
*/
	Milestone = 190;
	/*=============================================*/
	Pause();
	X = UfThold * V;
	Y = Radix * Radix;
	if (X*Y < One || X > Y) {
		if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
		else BadCond(Flaw, "");
			
# ifdef QP
		printf(" unbalanced range; UfThold * V = %.17Le\n\t%s\n",
			X, "is too far from 1.\n");
# else
		printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
			X, "is too far from 1.\n");
# endif
		}
	/*=============================================*/
	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))
# ifdef QP
			printf("  X / X  traps when X = %Lg\n", X);
# else
			printf("  X / X  traps when X = %g\n", X);
# endif
		else {
			V9 = (Y / X - Half) - Half;
			if (V9 == Zero) continue;
			if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
			else BadCond(Serious, "");
# ifdef QP
			printf("  X / X differs from 1 when X = %.17Le\n", X);
			printf("  instead, X / X - 1/2 - 1/2 = %.17Le .\n", V9);
# else
			printf("  X / X differs from 1 when X = %.17e\n", X);
			printf("  instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
# endif
			}
		}
	/*=============================================*/
	Milestone = 210;
	/*=============================================*/
	MyZero = Zero;
	printf("\n");
	printf("What message and/or values does Division by Zero produce?\n") ;
#ifndef NOPAUSE
	printf("This can interupt your program.  You can ");
	printf("skip this part if you wish.\n");
	printf("Do you wish to compute 1 / 0? ");
	fflush(stdout);
	read (KEYBOARD, ch, 8);
	if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif
		sigsave = Sigfpe;
		printf("    Trying to compute 1 / 0 produces ...");

#ifdef QP
		if (!setjmp(ovfl_buf)) printf("  %.7Le .\n", One / MyZero);
#else
		if (!setjmp(ovfl_buf)) printf("  %.7e .\n", One / MyZero);
#endif

#ifndef NOPAUSE
		}
	else printf("O.K.\n");
	printf("\nDo you wish to compute 0 / 0? ");
	fflush(stdout);
	read (KEYBOARD, ch, 80);
	if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif
		sigsave = Sigfpe;
		printf("\n    Trying to compute 0 / 0 produces ...");

#ifdef QP
		if (!setjmp(ovfl_buf)) printf("  %.7Le .\n", Zero / MyZero);
#else
		if (!setjmp(ovfl_buf)) printf("  %.7e .\n", Zero / MyZero);
#endif

#ifndef NOPAUSE
		}
	else printf("O.K.\n");
#endif
	/*=============================================*/
	Milestone = 220;
	/*=============================================*/
	Pause();
	printf("\n");
	{
		static char *msg[] = {
			"FAILUREs  encountered =",
			"SERIOUS DEFECTs  discovered =",
			"DEFECTs  discovered =",
			"FLAWs  discovered =" };
		int i;
		for(i = 0; i < 4; i++) if (ErrCnt[i])
			printf("The number of  %-29s %d.\n",
				msg[i], ErrCnt[i]);
		}
	printf("\n");
	if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
			+ ErrCnt[Flaw]) > 0) {
		if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
			Defect] == 0) && (ErrCnt[Flaw] > 0)) {
			printf("The arithmetic diagnosed seems ");
			printf("satisfactory though flawed.\n");
			}
		if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
			&& ( ErrCnt[Defect] > 0)) {
			printf("The arithmetic diagnosed may be acceptable\n");
			printf("despite inconvenient Defects.\n");
			}
		if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
			printf("The arithmetic diagnosed has ");
			printf("unacceptable serious defects.\n");
			}
		if (ErrCnt[Failure] > 0) {
			printf("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");
/* Modification for UCBTEST purposes - only perfect IEEE 754/854 counts. */
					(void) ucbpass( __FILE__ , __LINE__ );
					}
				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");
	(void) ucbfail( __FILE__ , __LINE__ );
	}

/*SPLIT subs.c
#include "paranoia.h"
*/

/* Sign */

GENERIC Sign (X)
GENERIC X;
{ return X >= 0. ? 1.0 : -1.0; }

/* Pause */

void
Pause()
{
	
#ifndef NOPAUSE
	char ch[8];
	printf("\nTo continue, press RETURN");
	fflush(stdout);
	read(KEYBOARD, ch, 8);
#endif
	printf("\nDiagnosis resumes after milestone Number %d", Milestone);
	printf("          Page: %d\n\n", PageNo);
	++Milestone;
	++PageNo;
	}

 /* TstCond */

void
TstCond (K, Valid, T)
int K, Valid;
char *T;
{ if (! Valid) { BadCond(K,T); printf(".\n"); } }

void
BadCond(K, T)
int K;
char *T;
{
	static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };

	ErrCnt [K] = ErrCnt [K] + 1;
	printf("%s:  %s", msg[K], T);
	}

/* Random */
/*  Random computes
     X = (Random1 + Random9)^5
     Random1 = X - FLOOR(X) + 0.000005 * X;
   and returns the new value of Random1
*/

GENERIC Random()
{
	GENERIC X, Y;
	
	X = Random1 + Random9;
	Y = X * X;
	Y = Y * Y;
	X = X * Y;
	Y = X - FLOOR(X);
	Random1 = Y + X * 0.000005;
	return(Random1);
	}

/* SqXMinX */

void
SqXMinX (ErrKind)
int ErrKind;
{
	GENERIC 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");
# ifdef QP
		printf("sqrt( %.17Le) - %.17Le  = %.17Le\n", X * X, X, OneUlp * SqEr);
# else
		printf("sqrt( %.17e) - %.17e  = %.17e\n", X * X, X, OneUlp * SqEr);
# endif
		printf("\tinstead of correct value 0 .\n");
		}
	}

/* NewD */

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;
	}

/* SR3750 */

void
SR3750()
{
	if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
		I = I + 1;
		X2 = SQRT(X * D);
		Y2 = (X2 - Z2) - (Y - Z2);
		X2 = X8 / (Y - Half);
		X2 = X2 - Half * X2 * X2;
		SqEr = (Y2 + Half) + (Half - X2);
		if (SqEr < MinSqEr) MinSqEr = SqEr;
		SqEr = Y2 - X2;
		if (SqEr > MaxSqEr) MaxSqEr = SqEr;
		}
	}

/* IsYeqX */

void
IsYeqX()
{
	if (Y != X) {
		if (N <= 0) {
			if (Z == Zero && Q <= Zero)
				printf("WARNING:  computing\n");
			else BadCond(Defect, "computing\n");
# ifdef QP
			printf("\t(%.17Le) ^ (%.17Le)\n", Z, Q);
			printf("\tyielded %.17Le;\n", Y);
			printf("\twhich compared unequal to correct %.17Le ;\n",
				X);
			printf("\t\tthey differ by %.17Le .\n", Y - X);
# else
			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);
# endif
			}
		N = N + 1; /* ... count discrepancies. */
		}
	}

/* SR3980 */

void
SR3980()
{
	do {
		Q = (GENERIC) I;
		Y = POW(Z, Q);
		IsYeqX();
		if (++I > M) break;
		X = Z * X;
		} while ( X < W );
	}

/* PrintIfNPositive */

void
PrintIfNPositive()
{
	if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
	}

/* TstPtUf */

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;
#ifdef QP
		printf("What the machine gets for (Z + Z) / Z is  %.17Le .\n",
			Q9);
#else
		printf("What the machine gets for (Z + Z) / Z is  %.17e .\n",
			Q9);
#endif
		if (ABS(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");
				}
			}
		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 = ");
# ifdef QP
			printf("%.17Le\n\tcompares different from  ", Z);
			if (Z != Random1) printf("Z * 1 = %.17Le ", Random1);
			if (! ((Z == Random2)
				|| (Random2 == Random1)))
				printf("1 * Z == %Lg\n", Random2);
			if (! (Z == V9)) printf("Z / 1 = %.17Le\n", V9);
# else
			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);
# endif
			if (Random2 != Random1) {
				ErrCnt [Defect] = ErrCnt [Defect] + 1;
				BadCond(Defect, "Multiplication does not commute!\n");
# ifdef QP
				printf("\tComparison alleges that 1 * Z = %.17Le\n",
					Random2);
				printf("\tdiffers from Z * 1 = %.17Le\n", Random1);
# else
				printf("\tComparison alleges that 1 * Z = %.17e\n",
					Random2);
				printf("\tdiffers from Z * 1 = %.17e\n", Random1);
# endif
				}
			Pause();
			}
		}
	}

void
notify(s)
char *s;
{
	printf("%s test appears to be inconsistent...\n", s);
	printf("   PLEASE NOTIFY KARPINKSI!\n");
	}

/*SPLIT msgs.c */

/* Instructions */

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);
	}

/* Heading */

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
	"\tPrecision:\tsingle;",
#else
	"\tPrecision:\tdouble;",
#endif
	"\tVersion:\t9 June 1986;",
	"\tComputer:\n",
	"\tCompiler:\n",
	"\tOptimization level:\n",
	"\tOther relevant compiler options:",
	0};

	msglist(head);
	}

/* Characteristics */

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);
	}
