Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 0 insertions, 2348 deletions

diff --git a/libm/ldouble/lparanoi.c b/libm/ldouble/lparanoi.c
deleted file mode 100644
index eb8fd25c7..000000000
--- a/libm/ldouble/lparanoi.c
+++ /dev/null
@@ -1,2348 +0,0 @@
-/*	A C version of Kahan's Floating Point Test "Paranoia"
-
-			Thos Sumner, UCSF, Feb. 1985
-			David Gay, BTL, Jan. 1986
-
-	This is a rewrite from the Pascal version by
-
-			B. A. Wichmann, 18 Jan. 1985
-
-	(and does NOT exhibit good C programming style).
-
-(C) Apr 19 1983 in BASIC version by:
-	Professor W. M. Kahan,
-	567 Evans Hall
-	Electrical Engineering & Computer Science Dept.
-	University of California
-	Berkeley, California 94720
-	USA
-
-converted to Pascal by:
-	B. A. Wichmann
-	National Physical Laboratory
-	Teddington Middx
-	TW11 OLW
-	UK
-
-converted to C by:
-
-	David M. Gay		and	Thos Sumner
-	AT&T Bell Labs			Computer Center, Rm. U-76
-	600 Mountainn Avenue		University of California
-	Murray Hill, NJ 07974		San Francisco, CA 94143
-	USA				USA
-
-with simultaneous corrections to the Pascal source (reflected
-in the Pascal source available over netlib).
-
-Reports of results on various systems from all the versions
-of Paranoia are being collected by Richard Karpinski at the
-same address as Thos Sumner.  This includes sample outputs,
-bug reports, and criticisms.
-
-You may copy this program freely if you acknowledge its source.
-Comments on the Pascal version to NPL, please.
-
-
-The C version catches signals from floating-point exceptions.
-If signal(SIGFPE,...) is unavailable in your environment, you may
-#define NOSIGNAL to comment out the invocations of signal.
-
-This source file is too big for some C compilers, but may be split
-into pieces.  Comments containing "SPLIT" suggest convenient places
-for this splitting.  At the end of these comments is an "ed script"
-(for the UNIX(tm) editor ed) that will do this splitting.
-
-By #defining Single when you compile this source, you may obtain
-a single-precision C version of Paranoia.
-
-
-The following is from the introductory commentary from Wichmann's work:
-
-The BASIC program of Kahan is written in Microsoft BASIC using many
-facilities which have no exact analogy in Pascal.  The Pascal
-version below cannot therefore be exactly the same.  Rather than be
-a minimal transcription of the BASIC program, the Pascal coding
-follows the conventional style of block-structured languages.  Hence
-the Pascal version could be useful in producing versions in other
-structured languages.
-
-Rather than use identifiers of minimal length (which therefore have
-little mnemonic significance), the Pascal version uses meaningful
-identifiers as follows [Note: A few changes have been made for C]:
-
-
-BASIC   C               BASIC   C               BASIC   C               
-
-   A                       J                       S    StickyBit
-   A1   AInverse           J0   NoErrors           T
-   B    Radix                    [Failure]         T0   Underflow
-   B1   BInverse           J1   NoErrors           T2   ThirtyTwo
-   B2   RadixD2                  [SeriousDefect]   T5   OneAndHalf
-   B9   BMinusU2           J2   NoErrors           T7   TwentySeven
-   C                             [Defect]          T8   TwoForty
-   C1   CInverse           J3   NoErrors           U    OneUlp
-   D                             [Flaw]            U0   UnderflowThreshold
-   D4   FourD              K    PageNo             U1
-   E0                      L    Milestone          U2
-   E1                      M                       V
-   E2   Exp2               N                       V0
-   E3                      N1                      V8
-   E5   MinSqEr            O    Zero               V9
-   E6   SqEr               O1   One                W
-   E7   MaxSqEr            O2   Two                X
-   E8                      O3   Three              X1
-   E9                      O4   Four               X8
-   F1   MinusOne           O5   Five               X9   Random1
-   F2   Half               O8   Eight              Y
-   F3   Third              O9   Nine               Y1
-   F6                      P    Precision          Y2
-   F9                      Q                       Y9   Random2
-   G1   GMult              Q8                      Z
-   G2   GDiv               Q9                      Z0   PseudoZero
-   G3   GAddSub            R                       Z1
-   H                       R1   RMult              Z2
-   H1   HInverse           R2   RDiv               Z9
-   I                       R3   RAddSub
-   IO   NoTrials           R4   RSqrt
-   I3   IEEE               R9   Random9
-
-   SqRWrng
-
-All the variables in BASIC are true variables and in consequence,
-the program is more difficult to follow since the "constants" must
-be determined (the glossary is very helpful).  The Pascal version
-uses Real constants, but checks are added to ensure that the values
-are correctly converted by the compiler.
-
-The major textual change to the Pascal version apart from the
-identifiersis that named procedures are used, inserting parameters
-wherehelpful.  New procedures are also introduced.  The
-correspondence is as follows:
-
-
-BASIC       Pascal
-lines 
-
-  90- 140   Pause
- 170- 250   Instructions
- 380- 460   Heading
- 480- 670   Characteristics
- 690- 870   History
-2940-2950   Random
-3710-3740   NewD
-4040-4080   DoesYequalX
-4090-4110   PrintIfNPositive
-4640-4850   TestPartialUnderflow
-
-=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
-
-Below is an "ed script" that splits para.c into 10 files
-of the form part[1-8].c, subs.c, and msgs.c, plus a header
-file, paranoia.h, that these files require.
-r paranoia.c
-$
-?SPLIT
-+,$w msgs.c
-.,$d
-?SPLIT
-.d
-+d
--,$w subs.c
--,$d
-?part8
-+d
-?include
-.,$w part8.c
-.,$d
--d
-?part7
-+d
-?include
-.,$w part7.c
-.,$d
--d
-?part6
-+d
-?include
-.,$w part6.c
-.,$d
--d
-?part5
-+d
-?include
-.,$w part5.c
-.,$d
--d
-?part4
-+d
-?include
-.,$w part4.c
-.,$d
--d
-?part3
-+d
-?include
-.,$w part3.c
-.,$d
--d
-?part2
-+d
-?include
-.,$w part2.c
-.,$d
-?SPLIT
-.d
-1,/^#include/-1d
-1,$w part1.c
-/Computed constants/,$d
-1,$s/^int/extern &/
-1,$s/^FLOAT/extern &/
-1,$s! = .*!;!
-/^Guard/,/^Round/s/^/extern /
-/^jmp_buf/s/^/extern /
-/^Sig_type/s/^/extern /
-a
-extern int sigfpe();
-.
-w paranoia.h
-q
-
-*/
-
-#include <stdio.h>
-#ifndef NOSIGNAL
-#include <signal.h>
-#endif
-#include <setjmp.h>
-
-#define Ldouble
-/*#define Single*/
-
-#ifdef Single
-#define NPRT 2
-extern double fabs(), floor(), log(), pow(), sqrt();
-#define FLOAT float
-#define FABS(x) (float)fabs((double)(x))
-#define FLOOR(x) (float)floor((double)(x))
-#define LOG(x) (float)log((double)(x))
-#define POW(x,y) (float)pow((double)(x),(double)(y))
-#define SQRT(x) (float)sqrt((double)(x))
-#define FSETUP sprec
-/*sprec() { }*/
-#else
-#ifdef Ldouble
-#define NPRT 6
-extern long double fabsl(), floorl(), logl(), powl(), sqrtl();
-#define FLOAT long double
-#define FABS(x) fabsl(x)
-#define FLOOR(x) floorl(x)
-#define LOG(x) logl(x)
-#define POW(x,y) powl(x,y)
-#define SQRT(x) sqrtl(x)
-#define FSETUP ldprec
-#else
-#define NPRT 4
-extern double fabs(), floor(), log(), pow(), sqrt();
-#define FLOAT double
-#define FABS(x) fabs(x)
-#define FLOOR(x) floor(x)
-#define LOG(x) log(x)
-#define POW(x,y) pow(x,y)
-#define SQRT(x) sqrt(x)
-/*double __sqrtdf2();
-#define SQRT(x) __sqrtdf2(x)
-*/
-#define FSETUP dprec
-/* dprec() { } */
-#endif
-#endif
-
-jmp_buf ovfl_buf;
-typedef int (*Sig_type)();
-Sig_type sigsave;
-
-#define KEYBOARD 0
-
-FLOAT Radix, BInvrse, RadixD2, BMinusU2;
-FLOAT Sign(), Random();
-
-/*Small floating point constants.*/
-FLOAT Zero = 0.0;
-FLOAT Half = 0.5;
-FLOAT One = 1.0;
-FLOAT Two = 2.0;
-FLOAT Three = 3.0;
-FLOAT Four = 4.0;
-FLOAT Five = 5.0;
-FLOAT Eight = 8.0;
-FLOAT Nine = 9.0;
-FLOAT TwentySeven = 27.0;
-FLOAT ThirtyTwo = 32.0;
-FLOAT TwoForty = 240.0;
-FLOAT MinusOne = -1.0;
-FLOAT OneAndHalf = 1.5;
-/*Integer constants*/
-int NoTrials = 20; /*Number of tests for commutativity. */
-#define False 0
-#define True 1
-
-/* Definitions for declared types 
-	Guard == (Yes, No);
-	Rounding == (Chopped, Rounded, Other);
-	Message == packed array [1..40] of char;
-	Class == (Flaw, Defect, Serious, Failure);
-	  */
-#define Yes 1
-#define No  0
-#define Chopped 2
-#define Rounded 1
-#define Other   0
-#define Flaw    3
-#define Defect  2
-#define Serious 1
-#define Failure 0
-typedef int Guard, Rounding, Class;
-typedef char Message;
-
-/* Declarations of Variables */
-int Indx;
-char ch[8];
-FLOAT AInvrse, A1;
-FLOAT C, CInvrse;
-FLOAT D, FourD;
-static FLOAT E0, E1, Exp2, E3, MinSqEr;
-FLOAT SqEr, MaxSqEr, E9;
-FLOAT Third;
-FLOAT F6, F9;
-FLOAT H, HInvrse;
-int I;
-FLOAT StickyBit, J;
-FLOAT MyZero;
-FLOAT Precision;
-FLOAT Q, Q9;
-FLOAT R, Random9;
-FLOAT T, Underflow, S;
-FLOAT OneUlp, UfThold, U1, U2;
-FLOAT V, V0, V9;
-FLOAT W;
-FLOAT X, X1, X2, X8, Random1;
-static FLOAT Y, Y1, Y2, Random2;
-FLOAT Z, PseudoZero, Z1, Z2, Z9;
-int ErrCnt[4];
-int fpecount;
-int Milestone;
-int PageNo;
-int M, N, N1;
-Guard GMult, GDiv, GAddSub;
-Rounding RMult, RDiv, RAddSub, RSqrt;
-int Break, Done, NotMonot, Monot, Anomaly, IEEE,
-		SqRWrng, UfNGrad;
-/* Computed constants. */
-/*U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
-/*U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
-
-/* floating point exception receiver */
-sigfpe()
-{
-	fpecount++;
-	printf("\n* * * FLOATING-POINT ERROR * * *\n");
-	fflush(stdout);
-	if (sigsave) {
-#ifndef NOSIGNAL
-		signal(SIGFPE, sigsave);
-#endif
-		sigsave = 0;
-		longjmp(ovfl_buf, 1);
-		}
-	abort();
-}
-
-
-FLOAT Ptemp;
-
-pnum( x )
-FLOAT *x;
-{
-char str[30];
-double d;
-unsigned short *p;
-int i;
-
-p = (unsigned short *)x;
-for( i=0; i<NPRT; i++ )
-	printf( "%04x ", *p++ & 0xffff );
-#ifdef Ldouble
-e64toasc( x, str, 20 );
-#else
-#ifdef Single
-e24toasc( x, str, 20 );
-#else
-e53toasc( x, str, 20 );
-#endif
-#endif
-printf( " = %s\n", str );
-/*
-d = *x;
-printf( " = %.16e\n", d );
-*/
-}
-
-
-
-main()
-{
-/* noexcept(); */
- FSETUP();
-	/* First two assignments use integer right-hand sides. */
-	Zero = 0;
-	One = 1;
-	Two = One + One;
-	Three = Two + One;
-	Four = Three + One;
-	Five = Four + One;
-	Eight = Four + Four;
-	Nine = Three * Three;
-	TwentySeven = Nine * Three;
-	ThirtyTwo = Four * Eight;
-	TwoForty = Four * Five * Three * Four;
-	MinusOne = -One;
-	Half = One / Two;
-	OneAndHalf = One + Half;
-	ErrCnt[Failure] = 0;
-	ErrCnt[Serious] = 0;
-	ErrCnt[Defect] = 0;
-	ErrCnt[Flaw] = 0;
-	PageNo = 1;
-	/*=============================================*/
-	Milestone = 0;
-	/*=============================================*/
-#ifndef NOSIGNAL
-	signal(SIGFPE, sigfpe);
-#endif
-	Instructions();
-	Pause();
-	Heading();
-	Pause();
-	Characteristics();
-	Pause();
-	History();
-	Pause();
-	/*=============================================*/
-	Milestone = 7;
-	/*=============================================*/
-	printf("Program is now RUNNING tests on small integers:\n");
-	
-	TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
-		   && (One > Zero) && (One + One == Two),
-			"0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
-	Z = - Zero;
-	if (Z == 0.0) {
-		U1 = 0.001;
-		Radix = 1;
-		TstPtUf();
-		}
-	else {
-		ErrCnt[Failure] = ErrCnt[Failure] + 1;
-		printf("Comparison alleges that -0.0 is Non-zero!\n");
-		}
-	TstCond (Failure, (Three == Two + One) && (Four == Three + One)
-		   && (Four + Two * (- Two) == Zero)
-		   && (Four - Three - One == Zero),
-		   "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
-	TstCond (Failure, (MinusOne == (0 - One))
-		   && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
-		   && (MinusOne + FABS(One) == Zero)
-		   && (MinusOne + MinusOne * MinusOne == Zero),
-		   "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
-	TstCond (Failure, Half + MinusOne + Half == Zero,
-		  "1/2 + (-1) + 1/2 != 0");
-	/*=============================================*/
-	/*SPLIT
-	part2();
-	part3();
-	part4();
-	part5();
-	part6();
-	part7();
-	part8();
-	}
-#include "paranoia.h"
-part2(){
-*/
-	Milestone = 10;
-	/*=============================================*/
-	TstCond (Failure, (Nine == Three * Three)
-		   && (TwentySeven == Nine * Three) && (Eight == Four + Four)
-		   && (ThirtyTwo == Eight * Four)
-		   && (ThirtyTwo - TwentySeven - Four - One == Zero),
-		   "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
-	TstCond (Failure, (Five == Four + One) &&
-			(TwoForty == Four * Five * Three * Four)
-		   && (TwoForty / Three - Four * Four * Five == Zero)
-		   && ( TwoForty / Four - Five * Three * Four == Zero)
-		   && ( TwoForty / Five - Four * Three * Four == Zero),
-		  "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
-	if (ErrCnt[Failure] == 0) {
-		printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
-		printf("\n");
-		}
-	printf("Searching for Radix and Precision.\n");
-	W = One;
-	do  {
-		W = W + W;
-		Y = W + One;
-		Z = Y - W;
-		Y = Z - One;
-		} while (MinusOne + FABS(Y) < Zero);
-	/*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
-	Precision = Zero;
-	Y = One;
-	do  {
-		Radix = W + Y;
-		Y = Y + Y;
-		Radix = Radix - W;
-		} while ( Radix == Zero);
-	if (Radix < Two) Radix = One;
-	printf("Radix = " );
-	pnum( &Radix );
-	if (Radix != 1) {
-		W = One;
-		do  {
-			Precision = Precision + One;
-			W = W * Radix;
-			Y = W + One;
-			} while ((Y - W) == One);
-		}
-	/*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
-			                              ...*/
-	U1 = One / W;
-	U2 = Radix * U1;
-	printf("Closest relative separation found is U1 = " );
-	pnum( &U1 );
-	printf("U2 = ");
-	pnum( &U2 );
-	printf("Recalculating radix and precision.");
-	
-	/*save old values*/
-	E0 = Radix;
-	E1 = U1;
-	E9 = U2;
-	E3 = Precision;
-	
-	X = Four / Three;
-	Third = X - One;
-	F6 = Half - Third;
-	X = F6 + F6;
-	X = FABS(X - Third);
-	if (X < U2) X = U2;
-	
-	/*... now X = (unknown no.) ulps of 1+...*/
-	do  {
-		U2 = X;
-		Y = Half * U2 + ThirtyTwo * U2 * U2;
-		Y = One + Y;
-		X = Y - One;
-		} while ( ! ((U2 <= X) || (X <= Zero)));
-	
-	/*... now U2 == 1 ulp of 1 + ... */
-	X = Two / Three;
-	F6 = X - Half;
-	Third = F6 + F6;
-	X = Third - Half;
-	X = FABS(X + F6);
-	if (X < U1) X = U1;
-	
-	/*... now  X == (unknown no.) ulps of 1 -... */
-	do  {
-		U1 = X;
-		Y = Half * U1 + ThirtyTwo * U1 * U1;
-		Y = Half - Y;
-		X = Half + Y;
-		Y = Half - X;
-		X = Half + Y;
-		} while ( ! ((U1 <= X) || (X <= Zero)));
-	/*... now U1 == 1 ulp of 1 - ... */
-	if (U1 == E1) printf("confirms closest relative separation U1 .\n");
-	else
-		{
-		printf("gets better closest relative separation U1 = " );
-		pnum( &U1 );
-		}
-	W = One / U1;
-	F9 = (Half - U1) + Half;
-	Radix = FLOOR(0.01 + U2 / U1);
-	if (Radix == E0) printf("Radix confirmed.\n");
-	else
-		{
-		printf("MYSTERY: recalculated Radix = " );
-		pnum( &Radix );
-		}
-	TstCond (Defect, Radix <= Eight + Eight,
-		   "Radix is too big: roundoff problems");
-	TstCond (Flaw, (Radix == Two) || (Radix == 10)
-		   || (Radix == One), "Radix is not as good as 2 or 10");
-	/*=============================================*/
-	Milestone = 20;
-	/*=============================================*/
-	TstCond (Failure, F9 - Half < Half,
-		   "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
-	X = F9;
-	I = 1;
-	Y = X - Half;
-	Z = Y - Half;
-	TstCond (Failure, (X != One)
-		   || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
-	X = One + U2;
-	I = 0;
-	/*=============================================*/
-	Milestone = 25;
-	/*=============================================*/
-	/*... BMinusU2 = nextafter(Radix, 0) */
-	BMinusU2 = Radix - One;
-	BMinusU2 = (BMinusU2 - U2) + One;
-	/* Purify Integers */
-	if (Radix != One)  {
-		X = - TwoForty * LOG(U1) / LOG(Radix);
-		Y = FLOOR(Half + X);
-		if (FABS(X - Y) * Four < One) X = Y;
-		Precision = X / TwoForty;
-		Y = FLOOR(Half + Precision);
-		if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
-		}
-	if ((Precision != FLOOR(Precision)) || (Radix == One)) {
-		printf("Precision cannot be characterized by an Integer number\n");
-		printf("of significant digits but, by itself, this is a minor flaw.\n");
-		}
-	if (Radix == One) 
-		printf("logarithmic encoding has precision characterized solely by U1.\n");
-	else
-		{
-		printf("The number of significant digits of the Radix is " );
-		pnum( &Precision );
-		}
-	TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
-		   "Precision worse than 5 decimal figures  ");
-	/*=============================================*/
-	Milestone = 30;
-	/*=============================================*/
-	/* Test for extra-precise subepressions */
-	X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
-	do  {
-		Z2 = X;
-		X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
-		} while ( ! ((Z2 <= X) || (X <= Zero)));
-	X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
-	do  {
-		Z1 = Z;
-		Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
-			+ One / Two)) + One / Two;
-		} while ( ! ((Z1 <= Z) || (Z <= Zero)));
-	do  {
-		do  {
-			Y1 = Y;
-			Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
-				)) + Half;
-			} while ( ! ((Y1 <= Y) || (Y <= Zero)));
-		X1 = X;
-		X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
-		} while ( ! ((X1 <= X) || (X <= Zero)));
-	if ((X1 != Y1) || (X1 != Z1)) {
-		BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
-		printf("respectively  " );
-		pnum( &X1 );
-		pnum( &Y1 );
-		pnum( &Z1 );
-		printf("are symptoms of inconsistencies introduced\n");
-		printf("by extra-precise evaluation of arithmetic subexpressions.\n");
-		notify("Possibly some part of this");
-		if ((X1 == U1) || (Y1 == U1) || (Z1 == U1))  printf(
-			"That feature is not tested further by this program.\n") ;
-		}
-	else  {
-		if ((Z1 != U1) || (Z2 != U2)) {
-			if ((Z1 >= U1) || (Z2 >= U2)) {
-				BadCond(Failure, "");
-				notify("Precision");
-				printf("\tU1 = " );
-				pnum( &U1 );
-				printf( "Z1 - U1 = " );
-				Ptemp = Z1-U1;
-				pnum( &Ptemp );
-				printf("\tU2 = " );
-				pnum( &U2 );
-				Ptemp = Z2-U2;
-				printf( "Z2 - U2 = " );
-				pnum( &Ptemp );
-				}
-			else {
-				if ((Z1 <= Zero) || (Z2 <= Zero)) {
-					printf("Because of unusual Radix = ");
-					pnum( &Radix );
-					printf(", or exact rational arithmetic a result\n");
-					printf("Z1 = " );
-					pnum( &Z1 );
-					printf( "or Z2 = " );
-					pnum( &Z2 );
-					notify("of an\nextra-precision");
-					}
-				if (Z1 != Z2 || Z1 > Zero) {
-					X = Z1 / U1;
-					Y = Z2 / U2;
-					if (Y > X) X = Y;
-					Q = - LOG(X);
-					printf("Some subexpressions appear to be calculated extra\n");
-					printf("precisely with about" );
-					Ptemp = Q / LOG(Radix);
-					pnum( &Ptemp );
-					printf( "extra B-digits, i.e.\n" );
-					Ptemp = Q / LOG(10.);
-					printf("roughly " );
-					pnum( &Ptemp );
-					printf( "extra significant decimals.\n");
-					}
-				printf("That feature is not tested further by this program.\n");
-				}
-			}
-		}
-	Pause();
-	/*=============================================*/
-	/*SPLIT
-	}
-#include "paranoia.h"
-part3(){
-*/
-	Milestone = 35;
-	/*=============================================*/
-	if (Radix >= Two) {
-		X = W / (Radix * Radix);
-		Y = X + One;
-		Z = Y - X;
-		T = Z + U2;
-		X = T - Z;
-		TstCond (Failure, X == U2,
-			"Subtraction is not normalized X=Y,X+Z != Y+Z!");
-		if (X == U2) printf(
-			"Subtraction appears to be normalized, as it should be.");
-		}
-	printf("\nChecking for guard digit in *, /, and -.\n");
-	Y = F9 * One;
-	Z = One * F9;
-	X = F9 - Half;
-	Y = (Y - Half) - X;
-	Z = (Z - Half) - X;
-	X = One + U2;
-	T = X * Radix;
-	R = Radix * X;
-	X = T - Radix;
-	X = X - Radix * U2;
-	T = R - Radix;
-	T = T - Radix * U2;
-	X = X * (Radix - One);
-	T = T * (Radix - One);
-	if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
-	else {
-		GMult = No;
-		TstCond (Serious, False,
-			"* lacks a Guard Digit, so 1*X != X");
-		}
-	Z = Radix * U2;
-	X = One + Z;
-	Y = FABS((X + Z) - X * X) - U2;
-	X = One - U2;
-	Z = FABS((X - U2) - X * X) - U1;
-	TstCond (Failure, (Y <= Zero)
-		   && (Z <= Zero), "* gets too many final digits wrong.\n");
-	Y = One - U2;
-	X = One + U2;
-	Z = One / Y;
-	Y = Z - X;
-	X = One / Three;
-	Z = Three / Nine;
-	X = X - Z;
-	T = Nine / TwentySeven;
-	Z = Z - T;
-	TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
-		"Division lacks a Guard Digit, so error can exceed 1 ulp\n\
-or  1/3  and  3/9  and  9/27 may disagree");
-	Y = F9 / One;
-	X = F9 - Half;
-	Y = (Y - Half) - X;
-	X = One + U2;
-	T = X / One;
-	X = T - X;
-	if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
-	else {
-		GDiv = No;
-		TstCond (Serious, False,
-			"Division lacks a Guard Digit, so X/1 != X");
-		}
-	X = One / (One + U2);
-	Y = X - Half - Half;
-	TstCond (Serious, Y < Zero,
-		   "Computed value of 1/1.000..1 >= 1");
-	X = One - U2;
-	Y = One + Radix * U2;
-	Z = X * Radix;
-	T = Y * Radix;
-	R = Z / Radix;
-	StickyBit = T / Radix;
-	X = R - X;
-	Y = StickyBit - Y;
-	TstCond (Failure, X == Zero && Y == Zero,
-			"* and/or / gets too many last digits wrong");
-	Y = One - U1;
-	X = One - F9;
-	Y = One - Y;
-	T = Radix - U2;
-	Z = Radix - BMinusU2;
-	T = Radix - T;
-	if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
-	else {
-		GAddSub = No;
-		TstCond (Serious, False,
-			"- lacks Guard Digit, so cancellation is obscured");
-		}
-	if (F9 != One && F9 - One >= Zero) {
-		BadCond(Serious, "comparison alleges  (1-U1) < 1  although\n");
-		printf("  subtration yields  (1-U1) - 1 = 0 , thereby vitiating\n");
-		printf("  such precautions against division by zero as\n");
-		printf("  ...  if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
-		}
-	if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
-		"     *, /, and - appear to have guard digits, as they should.\n");
-	/*=============================================*/
-	Milestone = 40;
-	/*=============================================*/
-	Pause();
-	printf("Checking rounding on multiply, divide and add/subtract.\n");
-	RMult = Other;
-	RDiv = Other;
-	RAddSub = Other;
-	RadixD2 = Radix / Two;
-	A1 = Two;
-	Done = False;
-	do  {
-		AInvrse = Radix;
-		do  {
-			X = AInvrse;
-			AInvrse = AInvrse / A1;
-			} while ( ! (FLOOR(AInvrse) != AInvrse));
-		Done = (X == One) || (A1 > Three);
-		if (! Done) A1 = Nine + One;
-		} while ( ! (Done));
-	if (X == One) A1 = Radix;
-	AInvrse = One / A1;
-	X = A1;
-	Y = AInvrse;
-	Done = False;
-	do  {
-		Z = X * Y - Half;
-		TstCond (Failure, Z == Half,
-			"X * (1/X) differs from 1");
-		Done = X == Radix;
-		X = Radix;
-		Y = One / X;
-		} while ( ! (Done));
-	Y2 = One + U2;
-	Y1 = One - U2;
-	X = OneAndHalf - U2;
-	Y = OneAndHalf + U2;
-	Z = (X - U2) * Y2;
-	T = Y * Y1;
-	Z = Z - X;
-	T = T - X;
-	X = X * Y2;
-	Y = (Y + U2) * Y1;
-	X = X - OneAndHalf;
-	Y = Y - OneAndHalf;
-	if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
-	  printf("Y2 = ");
-	  pnum( &Y2 );
-	  printf("Y1 = ");
-	  pnum( &Y1 );
-	  printf("U2 = ");
-	  pnum( &U2 );
-		X = (OneAndHalf + U2) * Y2;
-		Y = OneAndHalf - U2 - U2;
-		Z = OneAndHalf + U2 + U2;
-		T = (OneAndHalf - U2) * Y1;
-		X = X - (Z + U2);
-		StickyBit = Y * Y1;
-		S = Z * Y2;
-		T = T - Y;
-		Y = (U2 - Y) + StickyBit;
-		Z = S - (Z + U2 + U2);
-		StickyBit = (Y2 + U2) * Y1;
-		Y1 = Y2 * Y1;
-		StickyBit = StickyBit - Y2;
-		Y1 = Y1 - Half;
-		if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
-			&& ( StickyBit == Zero) && (Y1 == Half)) {
-			RMult = Rounded;
-			printf("Multiplication appears to round correctly.\n");
-			}
-		else	if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
-				&& (T < Zero) && (StickyBit + U2 == Zero)
-				&& (Y1 < Half)) {
-				RMult = Chopped;
-				printf("Multiplication appears to chop.\n");
-				}
-			else printf("* is neither chopped nor correctly rounded.\n");
-		if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
-		}
-	else printf("* is neither chopped nor correctly rounded.\n");
-	/*=============================================*/
-	Milestone = 45;
-	/*=============================================*/
-	Y2 = One + U2;
-	Y1 = One - U2;
-	Z = OneAndHalf + U2 + U2;
-	X = Z / Y2;
-	T = OneAndHalf - U2 - U2;
-	Y = (T - U2) / Y1;
-	Z = (Z + U2) / Y2;
-	X = X - OneAndHalf;
-	Y = Y - T;
-	T = T / Y1;
-	Z = Z - (OneAndHalf + U2);
-	T = (U2 - OneAndHalf) + T;
-	if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
-		X = OneAndHalf / Y2;
-		Y = OneAndHalf - U2;
-		Z = OneAndHalf + U2;
-		X = X - Y;
-		T = OneAndHalf / Y1;
-		Y = Y / Y1;
-		T = T - (Z + U2);
-		Y = Y - Z;
-		Z = Z / Y2;
-		Y1 = (Y2 + U2) / Y2;
-		Z = Z - OneAndHalf;
-		Y2 = Y1 - Y2;
-		Y1 = (F9 - U1) / F9;
-		if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
-			&& (Y2 == Zero) && (Y2 == Zero)
-			&& (Y1 - Half == F9 - Half )) {
-			RDiv = Rounded;
-			printf("Division appears to round correctly.\n");
-			if (GDiv == No) notify("Division");
-			}
-		else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
-			&& (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
-			RDiv = Chopped;
-			printf("Division appears to chop.\n");
-			}
-		}
-	if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
-	BInvrse = One / Radix;
-	TstCond (Failure, (BInvrse * Radix - Half == Half),
-		   "Radix * ( 1 / Radix ) differs from 1");
-	/*=============================================*/
-	/*SPLIT
-	}
-#include "paranoia.h"
-part4(){
-*/
-	Milestone = 50;
-	/*=============================================*/
-	TstCond (Failure, ((F9 + U1) - Half == Half)
-		   && ((BMinusU2 + U2 ) - One == Radix - One),
-		   "Incomplete carry-propagation in Addition");
-	X = One - U1 * U1;
-	Y = One + U2 * (One - U2);
-	Z = F9 - Half;
-	X = (X - Half) - Z;
-	Y = Y - One;
-	if ((X == Zero) && (Y == Zero)) {
-		RAddSub = Chopped;
-		printf("Add/Subtract appears to be chopped.\n");
-		}
-	if (GAddSub == Yes) {
-		X = (Half + U2) * U2;
-		Y = (Half - U2) * U2;
-		X = One + X;
-		Y = One + Y;
-		X = (One + U2) - X;
-		Y = One - Y;
-		if ((X == Zero) && (Y == Zero)) {
-			X = (Half + U2) * U1;
-			Y = (Half - U2) * U1;
-			X = One - X;
-			Y = One - Y;
-			X = F9 - X;
-			Y = One - Y;
-			if ((X == Zero) && (Y == Zero)) {
-				RAddSub = Rounded;
-				printf("Addition/Subtraction appears to round correctly.\n");
-				if (GAddSub == No) notify("Add/Subtract");
-				}
-			else printf("Addition/Subtraction neither rounds nor chops.\n");
-			}
-		else printf("Addition/Subtraction neither rounds nor chops.\n");
-		}
-	else printf("Addition/Subtraction neither rounds nor chops.\n");
-	S = One;
-	X = One + Half * (One + Half);
-	Y = (One + U2) * Half;
-	Z = X - Y;
-	T = Y - X;
-	StickyBit = Z + T;
-	if (StickyBit != Zero) {
-		S = Zero;
-		BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
-		}
-	StickyBit = Zero;
-	if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
-		&& (RMult == Rounded) && (RDiv == Rounded)
-		&& (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
-		printf("Checking for sticky bit.\n");
-		X = (Half + U1) * U2;
-		Y = Half * U2;
-		Z = One + Y;
-		T = One + X;
-		if ((Z - One <= Zero) && (T - One >= U2)) {
-			Z = T + Y;
-			Y = Z - X;
-			if ((Z - T >= U2) && (Y - T == Zero)) {
-				X = (Half + U1) * U1;
-				Y = Half * U1;
-				Z = One - Y;
-				T = One - X;
-				if ((Z - One == Zero) && (T - F9 == Zero)) {
-					Z = (Half - U1) * U1;
-					T = F9 - Z;
-					Q = F9 - Y;
-					if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
-						Z = (One + U2) * OneAndHalf;
-						T = (OneAndHalf + U2) - Z + U2;
-						X = One + Half / Radix;
-						Y = One + Radix * U2;
-						Z = X * Y;
-						if (T == Zero && X + Radix * U2 - Z == Zero) {
-							if (Radix != Two) {
-								X = Two + U2;
-								Y = X / Two;
-								if ((Y - One == Zero)) StickyBit = S;
-								}
-							else StickyBit = S;
-							}
-						}
-					}
-				}
-			}
-		}
-	if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
-	else printf("Sticky bit used incorrectly or not at all.\n");
-	TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
-			RMult == Other || RDiv == Other || RAddSub == Other),
-		"lack(s) of guard digits or failure(s) to correctly round or chop\n\
-(noted above) count as one flaw in the final tally below");
-	/*=============================================*/
-	Milestone = 60;
-	/*=============================================*/
-	printf("\n");
-	printf("Does Multiplication commute?  ");
-	printf("Testing on %d random pairs.\n", NoTrials);
-	Ptemp = 3.0;
-	Random9 = SQRT(Ptemp);
-	Random1 = Third;
-	I = 1;
-	do  {
-		X = Random();
-		Y = Random();
-		Z9 = Y * X;
-		Z = X * Y;
-		Z9 = Z - Z9;
-		I = I + 1;
-		} while ( ! ((I > NoTrials) || (Z9 != Zero)));
-	if (I == NoTrials) {
-		Random1 = One + Half / Three;
-		Random2 = (U2 + U1) + One;
-		Z = Random1 * Random2;
-		Y = Random2 * Random1;
-		Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
-			Three) * ((U2 + U1) + One);
-		}
-	if (! ((I == NoTrials) || (Z9 == Zero)))
-		BadCond(Defect, "X * Y == Y * X trial fails.\n");
-	else printf("     No failures found in %d integer pairs.\n", NoTrials);
-	/*=============================================*/
-	Milestone = 70;
-	/*=============================================*/
-	printf("\nRunning test of square root(x).\n");
-	TstCond (Failure, (Zero == SQRT(Zero))
-		   && (- Zero == SQRT(- Zero))
-		   && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
-	MinSqEr = Zero;
-	MaxSqEr = Zero;
-	J = Zero;
-	X = Radix;
-	OneUlp = U2;
-	SqXMinX (Serious);
-	X = BInvrse;
-	OneUlp = BInvrse * U1;
-	SqXMinX (Serious);
-	X = U1;
-	OneUlp = U1 * U1;
-	SqXMinX (Serious);
-	if (J != Zero) Pause();
-	printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
-	J = Zero;
-	X = Two;
-	Y = Radix;
-	if ((Radix != One)) do  {
-		X = Y;
-		Y = Radix * Y;
-		} while ( ! ((Y - X >= NoTrials)));
-	OneUlp = X * U2;
-	I = 1;
-	while (I < 10) {
-		X = X + One;
-		SqXMinX (Defect);
-		if (J > Zero) break;
-		I = I + 1;
-		}
-	printf("Test for sqrt monotonicity.\n");
-	I = - 1;
-	X = BMinusU2;
-	Y = Radix;
-	Z = Radix + Radix * U2;
-	NotMonot = False;
-	Monot = False;
-	while ( ! (NotMonot || Monot)) {
-		I = I + 1;
-		X = SQRT(X);
-		Q = SQRT(Y);
-		Z = SQRT(Z);
-		if ((X > Q) || (Q > Z)) NotMonot = True;
-		else {
-			Q = FLOOR(Q + Half);
-			if ((I > 0) || (Radix == Q * Q)) Monot = True;
-			else if (I > 0) {
-			if (I > 1) Monot = True;
-			else {
-				Y = Y * BInvrse;
-				X = Y - U1;
-				Z = Y + U1;
-				}
-			}
-			else {
-				Y = Q;
-				X = Y - U2;
-				Z = Y + U2;
-				}
-			}
-		}
-	if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
-	else {
-		BadCond(Defect, "");
-		printf("sqrt(X) is non-monotonic for X near " );
-		pnum( &Y );
-		}
-	/*=============================================*/
-	/*SPLIT
-	}
-#include "paranoia.h"
-part5(){
-*/
-	Milestone = 80;
-	/*=============================================*/
-	MinSqEr = MinSqEr + Half;
-	MaxSqEr = MaxSqEr - Half;
-	Y = (SQRT(One + U2) - One) / U2;
-	SqEr = (Y - One) + U2 / Eight;
-	if (SqEr > MaxSqEr) MaxSqEr = SqEr;
-	SqEr = Y + U2 / Eight;
-	if (SqEr < MinSqEr) MinSqEr = SqEr;
-	Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
-	SqEr = Y + U1 / Eight;
-	if (SqEr > MaxSqEr) MaxSqEr = SqEr;
-	SqEr = (Y + One) + U1 / Eight;
-	if (SqEr < MinSqEr) MinSqEr = SqEr;
-	OneUlp = U2;
-	X = OneUlp;
-	for( Indx = 1; Indx <= 3; ++Indx) {
-		Y = SQRT((X + U1 + X) + F9);
-		Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
-		Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
-		SqEr = (Y + Half) + Z;
-		if (SqEr < MinSqEr) MinSqEr = SqEr;
-		SqEr = (Y - Half) + Z;
-		if (SqEr > MaxSqEr) MaxSqEr = SqEr;
-		if (((Indx == 1) || (Indx == 3))) 
-			X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
-		else {
-			OneUlp = U1;
-			X = - OneUlp;
-			}
-		}
-	/*=============================================*/
-	Milestone = 85;
-	/*=============================================*/
-	SqRWrng = False;
-	Anomaly = False;
-	if (Radix != One) {
-		printf("Testing whether sqrt is rounded or chopped.\n");
-		D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
-	/* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
-		X = D / Radix;
-		Y = D / A1;
-		if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
-			Anomaly = True;
-			}
-		else {
-			X = Zero;
-			Z2 = X;
-			Y = One;
-			Y2 = Y;
-			Z1 = Radix - One;
-			FourD = Four * D;
-			do  {
-				if (Y2 > Z2) {
-					Q = Radix;
-					Y1 = Y;
-					do  {
-						X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
-						Q = Y1;
-						Y1 = X1;
-						} while ( ! (X1 <= Zero));
-					if (Q <= One) {
-						Z2 = Y2;
-						Z = Y;
-						}
-					}
-				Y = Y + Two;
-				X = X + Eight;
-				Y2 = Y2 + X;
-				if (Y2 >= FourD) Y2 = Y2 - FourD;
-				} while ( ! (Y >= D));
-			X8 = FourD - Z2;
-			Q = (X8 + Z * Z) / FourD;
-			X8 = X8 / Eight;
-			if (Q != FLOOR(Q)) Anomaly = True;
-			else {
-				Break = False;
-				do  {
-					X = Z1 * Z;
-					X = X - FLOOR(X / Radix) * Radix;
-					if (X == One) 
-						Break = True;
-					else
-						Z1 = Z1 - One;
-					} while ( ! (Break || (Z1 <= Zero)));
-				if ((Z1 <= Zero) && (! Break)) Anomaly = True;
-				else {
-					if (Z1 > RadixD2) Z1 = Z1 - Radix;
-					do  {
-						NewD();
-						} while ( ! (U2 * D >= F9));
-					if (D * Radix - D != W - D) Anomaly = True;
-					else {
-						Z2 = D;
-						I = 0;
-						Y = D + (One + Z) * Half;
-						X = D + Z + Q;
-						SR3750();
-						Y = D + (One - Z) * Half + D;
-						X = D - Z + D;
-						X = X + Q + X;
-						SR3750();
-						NewD();
-						if (D - Z2 != W - Z2) Anomaly = True;
-						else {
-							Y = (D - Z2) + (Z2 + (One - Z) * Half);
-							X = (D - Z2) + (Z2 - Z + Q);
-							SR3750();
-							Y = (One + Z) * Half;
-							X = Q;
-							SR3750();
-							if (I == 0) Anomaly = True;
-							}
-						}
-					}
-				}
-			}
-		if ((I == 0) || Anomaly) {
-			BadCond(Failure, "Anomalous arithmetic with Integer < ");
-			printf("Radix^Precision = " );
-			pnum( &W );
-			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 " ); 
-		Ptemp = MinSqEr - Half;
-		pnum( &Ptemp );
-		printf("to %.7e ulps.\n"); 
-		Ptemp = Half + MaxSqEr;
-		pnum( &Ptemp );
-		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 discrepancis 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 = FABS(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");
-			printf("PseudoZero that prints out as: " );
-			pnum( &PseudoZero );
-			X = - PseudoZero;
-			if (X <= Zero) {
-				printf("But -PseudoZero, which should be\n");
-				printf("positive, isn't; it prints out as " );
-				pnum( &X );
-				}
-			}
-		else {
-			BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
-			printf("value PseudoZero that prints out as ");
-			pnum( &PseudoZero );
-			}
-		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");
-			}
-		}
-	printf("Smallest strictly positive number found is E0 = ");
-	Pause();
-	pnum( &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 = ");
-			pnum( &UfThold );
-			printf(" coming down from " );
-			pnum( &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 = " );
-		pnum( &Q );
-		printf( "Y = " );
-		pnum( &Y );
-		printf ("|Q - Y| = " );
-		Ptemp = FABS(Q - Y2);
-		pnum( &Ptemp );
-		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");
-		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, "");
-			printf("X = " );
-			pnum( &X );
-			printf( "is not equal to Z = ");
-			pnum( &Z );
-			Z9 = X - Z;
-			printf("yet X - Z yields " );
-			pnum( &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");
-			printf("X / Z = 1 + ");
-			Ptemp = (X / Z - Half) - Half;
-			pnum( &Ptemp );
-			}
-		}
-	printf("The Underflow threshold is ");
-	pnum( &UfThold );
-	printf("below which 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 - One;
-	printf("Since underflow occurs below the threshold\n");
-	printf("UfThold = ");
-	pnum( &HInvrse );
-	printf( ") ^ (Y=" );
-	pnum( &Y );
-	printf( ")\nonly underflow " );
-	printf("should afflict the expression HInvrse^(Y+1).\n");
-	pnum( &HInvrse );
-	pnum( &Y2 );
-	V9 = POW(HInvrse, Y2);
-	printf("actually calculating yields: ");
-	pnum( &V9 );
-	if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
-		BadCond(Serious, "this is not between 0 and underflow\n");
-		printf("   threshold = ");
-		pnum( &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 = ");
-		pnum( &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) = ");
-	pnum( &Exp2 );
-	printf( "as X -> 1.\n");
-	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");
-			Ptemp = POW(X,Z);
-			pnum(&Ptemp);
-			printf("for (1 + (" );
-			pnum( &V9 );
-			printf( ") ^ (" );
-			pnum( &Z );
-			printf(") differs from correct value by ");
-			pnum( &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");
-	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");
-	printf("Trying it on Y = " );
-	pnum( &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 0
-/*  this doesn't handle infinity. */
-	if (Z != Y) {
-		BadCond(Serious, "");
-		printf("overflow past " );
-		pnum( &Y );
-		printf( "shrinks to " );
-		pnum( &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;
-	printf("Overflow threshold is V  = " );
-	pnum( &V );
-	if (I)
-		{
-		printf("Overflow saturates at V0 = " );
-		pnum( &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 = " );
-	pnum( &V9 );
-	V9 = V / One;
-	printf("                           nor for V / 1 = " );
-	pnum( &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("+-" );
-		pnum( &V );
-		printf( ", +- " );
-		pnum( &V0 );
-		printf( "and +- " );
-		pnum( &UfThold );
-		printf( "are confused by Overflow." );
-		}
-	/*=============================================*/
-	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, "");
-				printf("Comparison alleges that what prints as Z =" );
-				pnum( &Z );
-				printf(" is too far from sqrt(Z) ^ 2 = ");
-				pnum( &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 < W) BadCond(Serious, "");
-			else BadCond(Defect, "");
-			printf("Comparison alleges that Z = ");
-			pnum( &Z );
-			printf(" is too far from sqrt(Z) ^ 2 " );
-			pnum( &Y );
-			}
-		}
-	/*=============================================*/
-	/*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, "");
-			
-		printf(" unbalanced range; UfThold * V = " );
-		pnum( &X );
-		printf( "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 = ");
-			pnum( &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 =");
-			pnum( &X );
-			printf("  instead, X / X - 1/2 - 1/2 = ");
-			pnum( &V9 );
-			}
-		}
-	/*=============================================*/
-	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 ...");
-		if (!setjmp(ovfl_buf))
-			{
-			Ptemp = One / MyZero;
-			pnum( &Ptemp );
-			}
-#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 ...");
-		if (!setjmp(ovfl_buf))
-			{
-			Ptemp = Zero / MyZero;
-			pnum( &Ptemp );
-			}
-#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");
-				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");
-	}
-
-/*SPLIT subs.c
-#include "paranoia.h"
-*/
-
-/* Sign */
-
-FLOAT Sign (X)
-FLOAT X;
-{ return X >= 0. ? 1.0 : -1.0; }
-
-/* Pause */
-
-Pause()
-{
-	char ch[8];
-	
-#ifndef NOPAUSE
-	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 */
-
-TstCond (K, Valid, T)
-int K, Valid;
-char *T;
-{ if (! Valid) { BadCond(K,T); printf(".\n"); } }
-
-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
-*/
-
-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);
-	}
-
-/* SqXMinX */
-
-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( ");
-		Ptemp = X * X;
-		pnum( &Ptemp );
-		printf( ") - " );
-		pnum( &X );
-		printf("  =  " );
-		Ptemp = OneUlp * SqEr;
-		pnum( &Ptemp );
-		printf("\tinstead of correct value 0 .\n");
-		}
-	}
-
-/* NewD */
-
-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 */
-
-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 */
-
-IsYeqX()
-{
-	if (Y != X) {
-		if (N <= 0) {
-			if (Z == Zero && Q <= Zero)
-				printf("WARNING:  computing\n");
-			else BadCond(Defect, "computing\n");
-			printf("\t(");
-			pnum( &Z );
-			printf( ") ^ (" );
-			pnum( &Q );
-			printf("\tyielded " );
-			pnum( &Y );
-			printf("\twhich compared unequal to correct " );
-			pnum( &X );
-			printf("\t\tthey differ by " );
-			Ptemp = Y - X;
-			pnum( &Ptemp );
-			}
-		N = N + 1; /* ... count discrepancies. */
-		}
-	}
-
-/* SR3980 */
-
-SR3980()
-{
-	do {
-		Q = (FLOAT) I;
-		Y = POW(Z, Q);
-		IsYeqX();
-		if (++I > M) break;
-		X = Z * X;
-		} while ( X < W );
-	}
-
-/* PrintIfNPositive */
-
-PrintIfNPositive()
-{
-	if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
-	}
-
-/* TstPtUf */
-
-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 " );
-		pnum( &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");
-				}
-			}
-		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 = ");
-			pnum( &Z );
-			printf("\tcompares different from  ");
-			if (Z != Random1)
-				{
-				printf("Z * 1 = " );
-				pnum( &Random1 );
-				}
-			if (! ((Z == Random2)
-				|| (Random2 == Random1)))
-				{
-				printf("1 * Z == " );
-				pnum( &Random2 );
-				}
-			if (! (Z == V9))
-				{
-				printf("Z / 1 = ");
-				pnum( &V9 );
-				}
-			if (Random2 != Random1) {
-				ErrCnt [Defect] = ErrCnt [Defect] + 1;
-				BadCond(Defect, "Multiplication does not commute!\n");
-				printf("\tComparison alleges that 1 * Z = ");
-				pnum( &Random2 );
-				printf("\tdiffers from Z * 1 = ");
-				pnum( &Random1 );
-				}
-			Pause();
-			}
-		}
-	}
-
-notify(s)
-char *s;
-{
-	printf("%s test appears to be inconsistent...\n", s);
-	printf("   PLEASE NOTIFY KARPINKSI!\n");
-	}
-
-/*SPLIT msgs.c */
-
-/* Instructions */
-
-msglist(s)
-char **s;
-{ while(*s) printf("%s\n", *s++); }
-
-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 */
-
-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:\t27 January 1986;",
-	"\tComputer:\n",
-	"\tCompiler:\n",
-	"\tOptimization level:\n",
-	"\tOther relevant compiler options:",
-	0};
-
-	msglist(head);
-	}
-
-/* Characteristics */
-
-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);
-	}
-
-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);
-	}