if (Diff(Sum(a, basem1), a) != 0.0) {		if (f_radix == 2) basem1=0.375;		else basem1=1.0;		if (Diff(Sum(a, basem1), a) != 0.0) irnd=2; /* away from 0 */		else irnd=1; /* to nearest */	} else irnd=0; /* towards 0 */	basem1=Diff(base, 0.5);	if (Diff(Diff(-a, basem1), -a) != 0.0) {		if (f_radix == 2) basem1=0.375;		else basem1=1.0;		if (Diff(Diff(-a, basem1), -a) != 0.0) mrnd=2; /* away from 0*/		else mrnd=1; /* to nearest */	} else mrnd=0; /* towards 0 */	f_rounds=4; /* Unknown rounding */	if (irnd==0 && mrnd==0) f_rounds=0; /* zero = chops */	if (irnd==1 && mrnd==1) f_rounds=1; /* nearest */	if (irnd==2 && mrnd==0) f_rounds=2; /* +inf */	if (irnd==0 && mrnd==2) f_rounds=3; /* -inf */	if (f_rounds != 4) {		Vprintf("%sArithmetic rounds towards ", co);		switch (f_rounds) {		      case 0: Vprintf("zero (i.e. it chops)"); break;		      case 1: Vprintf("nearest"); break;		      case 2: Vprintf("+infinity"); break;		      case 3: Vprintf("-infinity"); break;		      default: Vprintf("???"); break;		}		Vprintf("%s\n", oc);	} else { /* Hmm, try to give some help here: */		Vprintf("%sArithmetic rounds oddly: %s\n", co, oc);		Vprintf("%s    Negative numbers %s%s\n",			co, mrnd==0 ? "towards zero" :			    mrnd==1 ? "to nearest" :				      "away from zero",			oc);		Vprintf("%s    Positive numbers %s%s\n",			co, irnd==0 ? "towards zero" :			    irnd==1 ? "to nearest" :				      "away from zero",			oc);	}	/* An extra goody */	if (f_radix == 2 && f_rounds == 1) {		if (Diff(Sum(a, 1.0), a) != 0.0) {			Vprintf("%s   Tie breaking rounds up%s\n", co, oc);		} else if (Diff(Sum(a, 3.0), a) == 4.0) {			Vprintf("%s   Tie breaking rounds to even%s\n", co, oc);		} else {			Vprintf("%s   Tie breaking rounds down%s\n", co, oc);		}	}#ifdef PASS1 /* only for FLT */	if (F) i_define("FLT", "_ROUNDS", (long) f_rounds, (long) F_ROUNDS);#endif	/* Various flavours of epsilon ************************************/	negeps=f_mant_dig+f_mant_dig;	basein=1.0/base;	a=1.0;	for(i=1; i<=negeps; i++) a*=basein;	b=a;	while (Diff(Diff(1.0, a), 1.0) == 0.0) {		a*=base;		negeps--;	}	negeps= -negeps;	Vprintf("%sSmallest x such that 1.0-base**x != 1.0 = %d%s\n",		co, negeps, oc);	epsneg=a;	if ((f_radix!=2) && irnd) {	/*	a=(a*(1.0+a))/(1.0+1.0); => */		a=Div(Mul(a, Sum(1.0, a)), Sum(1.0, 1.0));	/*	if ((1.0-a)-1.0 != 0.0) epsneg=a; => */		if (Diff(Diff(1.0, a), 1.0) != 0.0) epsneg=a;	}	Vprintf("%sSmall x such that 1.0-x != 1.0 = %s%s\n",		co, f_rep(digs, (Long_double) epsneg), oc);	/* it may not be the smallest */	if (V) F_check(digs, (Long_double) epsneg);	Unexpected(19);	machep= -f_mant_dig-f_mant_dig;	a=b;	while (Diff(Sum(1.0, a), 1.0) == 0.0) { a*=base; machep++; }	Vprintf("%sSmallest x such that 1.0+base**x != 1.0 = %d%s\n",		co, machep, oc);	f_epsilon=a;	if ((f_radix!=2) && irnd) {	/*	a=(a*(1.0+a))/(1.0+1.0); => */		a=Div(Mul(a, Sum(1.0, a)), Sum(1.0, 1.0));	/*	if ((1.0+a)-1.0 != 0.0) f_epsilon=a; => */		if (Diff(Sum(1.0, a), 1.0) != 0.0) f_epsilon=a;	}	Vprintf("%sSmallest x such that 1.0+x != 1.0 = %s%s\n",		co, f_rep(digs, (Long_double) f_epsilon), oc);	/* Possible loss of precision warnings here from non-stdc compilers: */	if (F) f_define(Fname, "_EPSILON", digs, (Long_double) f_epsilon, MARK);	if (V || F) F_check(digs, (Long_double) f_epsilon);	Unexpected(20);	if (F) Validate(digs, (Long_double) f_epsilon, (Long_double) F_EPSILON,			f_epsilon == Self(F_EPSILON));	Unexpected(21);	/* Extra chop info *************************************************/	if (f_rounds == 0) {		if (Diff(Mul(Sum(1.0,f_epsilon),1.0),1.0) !=  0.0) {			Vprintf("%sAlthough arithmetic chops, it uses guard digits%s\n", co, oc);		}	}	/* Size of and minimum normalised exponent ************************/	y=0; i=0; k=1; z=basein; z1=(1.0+f_epsilon)/base;	/* Coarse search for the largest power of two */	if (setjmp(lab)==0) { /* for underflow trap */ /* Yields i, k, y, y1 */		do {			y=z; y1=z1;			z=Mul(y,y); z1=Mul(z1, y);			a=Mul(z,1.0);			z2=Div(z1,y);			if (z2 != y1) break;			if ((Sum(a,a) == 0.0) || (fabs(z) >= y)) break;			i++;			k+=k;		} while(1);	} else {		Vprintf("%s%s underflow generates a trap%s\n", co, Thing, oc);	}	Unexpected(22);	if (f_radix != 10) {		iexp=i+1; /* for the sign */		mx=k+k;	} else {		iexp=2;		iz=f_radix;		while (k >= iz) { iz*=f_radix; iexp++; }		mx=iz+iz-1;	}	/* Fine tune starting with y and y1 */	if (setjmp(lab)==0) { /* for underflow trap */ /* Yields k, f_min */		do {			f_min=y; z1=y1;			y=Div(y,base); y1=Div(y1,base);			a=Mul(y,1.0);			z2=Mul(y1,base);			if (z2 != z1) break;			if ((Sum(a,a) == 0.0) || (fabs(y) >= f_min)) break;			k++;		} while (1);	}	Unexpected(23);	f_min_exp=(-k)+1;	if ((mx <= k+k-3) && (f_radix != 10)) { mx+=mx; iexp+=1; }	Vprintf("%sNumber of bits used for exponent = %d%s\n", co, iexp, oc);	Vprintf("%sMinimum normalised exponent = %d%s\n", co, f_min_exp, oc);	if (F) i_define(Fname, "_MIN_EXP", (long) f_min_exp, (long) F_MIN_EXP);	if (setjmp(lab)==0) {		Vprintf("%sMinimum normalised positive number = %s%s\n",			co, f_rep(digs, (Long_double) f_min), oc);	} else {		eek_a_bug("printf can't print the smallest normalised number");		printf("\n");	}	Unexpected(24);	/* Possible loss of precision warnings here from non-stdc compilers: */	if (setjmp(lab) == 0) {		if (F) f_define(Fname, "_MIN", digs, (Long_double) f_min, MARK);		if (V || F) F_check(digs, (Long_double) f_min);	} else {		eek_a_bug("xxx_MIN caused a trap");		printf("\n");	}	if (setjmp(lab) == 0) {		if (F) Validate(digs, (Long_double) f_min, (Long_double) F_MIN,				f_min == Self(F_MIN));	} else {		printf("%s*** Verify failed for above #define!\n    %s %s\n\n",		       co, "Compiler has an unusable number for value", oc);		bugs++;	}	Unexpected(25);	a=1.0; f_min_10_exp=0;	while (a > f_min*10.0) { a/=10.0; f_min_10_exp--; }	if (F) i_define(Fname, "_MIN_10_EXP", (long) f_min_10_exp,			(long) F_MIN_10_EXP);	/* Minimum exponent ************************************************/	if (setjmp(lab)==0) { /* for underflow trap */ /* Yields xminner */		do {			xminner=y;			y=Div(y,base);			a=Mul(y,1.0);			if ((Sum(a,a) == 0.0) || (fabs(y) >= xminner)) break;		} while (1);	}	Unexpected(26);	if (xminner != 0.0 && xminner != f_min) {		normal= 0;		Vprintf("%sThe smallest numbers are not kept normalised%s\n",			co, oc);		if (setjmp(lab)==0) {		    Vprintf("%sSmallest unnormalised positive number = %s%s\n",			    co, f_rep(digs, (Long_double) xminner), oc);		    if (V) F_check(digs, (Long_double) xminner);		} else {			eek_a_bug("printf can't print the smallest unnormalised number.");			printf("\n");		}		Unexpected(27);	} else {		normal= 1;		Vprintf("%sThe smallest numbers are normalised%s\n", co, oc);	}	/* Maximum exponent ************************************************/	f_max_exp=2; f_max=1.0; newxmax=base+1.0;	inf=0; trap=0;	while (f_max<newxmax) {		f_max=newxmax;		if (setjmp(lab) == 0) { /* Yields inf, f_max_exp */			newxmax=Mul(newxmax, base);		} else {			trap=1;			break;		}		if (Div(newxmax, base) != f_max) {			inf=1; /* ieee infinity */			break;		}		f_max_exp++;	}	Unexpected(28);	if (trap) {		Vprintf("%s%s overflow generates a trap%s\n", co, Thing, oc);	}	if (inf) Vprintf("%sThere is an 'infinite' value%s\n", co, oc);	Vprintf("%sMaximum exponent = %d%s\n", co, f_max_exp, oc);	if (F) i_define(Fname, "_MAX_EXP", (long) f_max_exp, (long) F_MAX_EXP);	/* Largest number ***************************************************/	f_max=Diff(1.0, epsneg);	if (Mul(f_max,1.0) != f_max) f_max=Diff(1.0, Mul(base,epsneg));	for (i=1; i<=f_max_exp; i++) f_max=Mul(f_max, base);	if (setjmp(lab)==0) {		Vprintf("%sMaximum number = %s%s\n",			co, f_rep(digs, (Long_double) f_max), oc);	} else {		eek_a_bug("printf can't print the largest double.");		printf("\n");	}	if (setjmp(lab)==0) {	/* Possible loss of precision warnings here from non-stdc compilers: */		if (F) f_define(Fname, "_MAX", digs, (Long_double) f_max, MARK);		if (V || F) F_check(digs, (Long_double) f_max);	} else {		eek_a_bug("xxx_MAX caused a trap");		printf("\n");	}	if (setjmp(lab)==0) {		if (F) Validate(digs, (Long_double) f_max, (Long_double) F_MAX,				f_max == Self(F_MAX));	} else {		printf("%s*** Verify failed for above #define!\n    %s %s\n\n",		       co, "Compiler has an unusable number for value", oc);		bugs++;	}	Unexpected(29);	a=1.0; f_max_10_exp=0;	while (a < f_max/10.0) { a*=10.0; f_max_10_exp++; }	if (F) i_define(Fname, "_MAX_10_EXP", (long) f_max_10_exp,			(long) F_MAX_10_EXP);	/* Hidden bit + sanity check ****************************************/	if (f_radix != 10) {		hidden=0;		mantbits=floor_log(2, (Long_double)f_radix)*f_mant_dig;		if (mantbits+iexp == (int)sizeof(Number)*bits_per_byte) {			hidden=1;			Vprintf("%sArithmetic uses a hidden bit%s\n", co, oc);		} else if (mantbits+iexp+1 == (int)sizeof(Number)*bits_per_byte) {			Vprintf("%sArithmetic doesn't use a hidden bit%s\n",				co, oc);		} else {			printf("\n%s%s\n    %s %s %s!%s\n\n",			       co,			       "*** Something fishy here!",			       "Exponent size + mantissa size doesn't match",			       "with the size of a", thing,			       oc);		}		if (hidden && f_radix == 2 && f_max_exp+f_min_exp==3) {			Vprintf("%sIt looks like %s length IEEE format%s\n",				co, f_mant_dig==24 ? "single" :				    f_mant_dig==53 ? "double" :				    f_mant_dig >53 ? "extended" :						"some", oc);			if (f_rounds != 1 || normal) {				Vprintf("%s   though ", co);				if (f_rounds != 1) {					Vprintf("the rounding is unusual");					if (normal) Vprintf(" and ");				}				if (normal) Vprintf("the normalisation is unusual");				Vprintf("%s\n", oc);			}		} else {			Vprintf("%sIt doesn't look like IEEE format%s\n",				co, oc);		}	}	printf("\n"); /* regardless of verbosity */	return f_mant_dig;}Procedure EPROP(fprec, dprec, lprec) int fprec, dprec, lprec; {	/* See if expressions are evaluated in extended precision.	   Some compilers optimise even if you don't want it,	   and then this function fails to produce the right result.	   We try to diagnose this if it happens.	*/	volatile int eprec;	volatile double a, b, base, old;	volatile Number d, oldd, dbase, one, zero;	volatile int bad=0;	/* Size of mantissa **************************************/	a=1.0;	if (setjmp(lab) == 0) { /* Yields nothing */		do { old=a; a=a+a; }		while ((((a+1.0)-a)-1.0) == 0.0 && a>old);	} else bad=1;	/* Avoid the comparison if bad is set,	   to avoid trouble on the convex.  */	if (!bad && (a <= old)) bad=1;	if (!bad) {		b=1.0;		if (setjmp(lab) == 0) { /* Yields nothing */			do { old=b; b=b+b; }			while ((base=((a+b)-a)) == 0.0 && b>old);			if (b <= old) bad=1;		} else bad=1;	}	if (!bad) {		eprec=0; d=1.0; dbase=base; one=1.0; zero=0.0;		if (setjmp(lab) == 0) { /* Yields nothing */			do { eprec++; oldd=d; d=d*dbase; }			while ((((d+one)-d)-one) == zero && d>oldd);			if (d <= oldd) bad=1;		} else bad=1;	}	Unexpected(30);	if (bad) {	  Vprintf("%sCan't determine precision for %s expressions:\n%s%s\n", 		 co, thing, "   check that you compiled without optimisation!",		 oc);	} else if (eprec==dprec) {	  Vprintf("%s%s expressions are evaluated in double precision%s\n",		  co, Thing, oc);	} else if (eprec==fprec) {	  Vprintf("%s%s expressions are evaluated in float precision%s\n",		  co, Thing, oc);	} else if (eprec==lprec) {	  Vprintf("%s%s expressions are evaluated in long double precision%s\n",		  co, Thing, oc);	} else {		Vprintf("%s%s expressions are evaluated in a %s %s %d %s%s\n",			co, Thing, eprec>dprec ? "higher" : "lower",			"precision than double,\n   using",			eprec, "base digits",		        oc);	}}#else /* Number */#ifdef FPROP/* ARGSUSED */int FPROP(bits_per_byte) int bits_per_byte; {	return 0;}#endif#ifdef EPROP/* ARGSUSED */Procedure EPROP(fprec, dprec, lprec) int fprec, dprec, lprec; {}#endif#endif /* ifdef Number */#ifdef PASS3#undef PASS#endif#ifdef PASS2#undef PASS2#define PASS3 1#endif#ifdef PASS1#undef PASS1#define PASS2 1#endif/* If your C compiler doesn't accept the next #include,   replace __FILE__ with the file name - and get a new C compiler... */#ifdef PASS#include "hard-params.c"#endif