⟦cedff05a8⟧

TextFile

/* Determine some properties of C types on your machine/compiler
   Author: Steven Pemberton, CWI, Amsterdam; steven@cwi.nl
   Bugfixes and upgrades gratefully received.

   The program only works if overflows are ignored by the C system or are
   catchable by signal().

   If your C system is not unix but does have signal/setjmp, compile with
	cc -DSIGNAL hard.c
   otherwise with
	cc hard.c
   Don't use any optimisation flags: the program won't work if you do.
   Some compilers need a -f flag for floating point.

   You may need to add some calls to signal() for other sorts of exception
   on your machine than SIGFPE, and SIGOVER. See lines beginning #ifdef
   SIGNAL later in the program.

   Output is produced as C style comments so that the program can be used to
   produce a .h file with minimum upheaval.

   I apologise unreservedly for the contorted use of the preprocessor...

   If your C preprocessor doesn't have the predefined __FILE__ macro, and
   you want to call this file anything other than hard.c, change the
   following #define command accordingly.
*/


#ifndef __FILE__
#define __FILE__ "hard.c"
#endif

#ifndef PASS
#define PASS 1

#include <stdio.h>

#ifdef SysV
#include	<sys/utsname.h>
#include	<time.h>
#endif

#ifndef SIGNAL
#ifdef unix
#define SIGNAL
#endif /*unix*/
#endif /*SIGNAL*/

#ifdef SIGNAL

#include "../install/signal.h"

#include <setjmp.h>

	jmp_buf lab;
	overflow(sig) int sig; { /* what to do on overflow/underflow */
		(void) signal(sig, overflow);
		longjmp(lab, 1);
	}

#else /*!SIGNAL*/
	/* Dummy routines instead */
	int lab=1;
	int setjmp(lab) int lab; { return(0); }

#endif /*SIGNAL*/

main() {
	int bits; /* the number of bits per unit returned by sizeof() */
	int dprec, eprec, basic(), fprop(), dprop(), eprop();
	char *malloc();
	unsigned int size;
	long total;
	void fault();

#ifdef SysV
	struct utsname name;
	long	tloc;
#endif

#ifdef SIGNAL
#ifdef SIGFPE
	(void) signal(SIGFPE, overflow);
#endif
#ifdef SIGOVER
	(void) signal(SIGOVER, overflow);
#endif
/* Add more calls as necessary */
#endif /* SIGNAL */

	if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }

#ifdef SysV
	/* show the machine name and date at the top of the output for reference */
	printf("/\* Machine: ");
	uname(&name);
	printf("%s %s %s %s %s *\/\n",
		name.sysname, 
		name.nodename, 
		name.release, 
		name.version, 
		name.machine);
	tloc = time( (long *) 0);
	/* precision for date string required to avoid \n returned from asctime() */
	printf("/\* Date: %.24s *\/\n", asctime(localtime(&tloc)) );
#endif

	bits=  basic();
	(void) fprop(bits);
	dprec= dprop(bits);
	eprec= eprop();
	printf("\n");
	if (eprec!=dprec)
		printf("/\* Expressions are evaluated in a %s %s %d %s *\/\n",
			eprec>dprec ? "higher" : "lower (tut!)",
			"precision than double, using", eprec, "base digits");
	else printf("/\* Expressions are evaluated in double precision *\/\n");

	/* An extra goody: the approximate amount of data-space */
	/* Allocate store until no more available */
	size=1<<((bits*sizeof(int))-2);
	/* 
	 * NEW VERSION of mallocatable
	 *
 	 * fix to avoid troubles with some implementations of malloc(),
 	 * we use signal(SIGSEGV) with setjmp() and longjmp().
 	 *
	 * See also the above function fault().
	 *
 	 * (Philippe Dax : May 1988
 	 */
	for (signal(SIGSEGV, fault), total=0; size; size >>=1) {
		if (!setjmp(lab)) {
			for (;;) {
				if (malloc(size) == NULL)
					break;
				total += (size/2);
			}
		}
	}
	printf("Espace memoire disponible pour l'utilisateur = %ld Koctets\n", (total+511)/512);
}

/*
 * fix to avoid troubles with some implementations of malloc(),
 * we use signal(SIGSEGV) with setjmp() and longjmp().
 *
 * (Philippe Dax : May 1988
 */
void fault()
{
	static int val = 0;

	val++;
	signal(SIGSEGV, fault);
	longjmp(lab, val);
}

int basic() {
	/* The properties of the basic types.
	   Returns number of bits per sizeof unit */

	char c; int bits;

	if (setjmp(lab)!=0) { printf("\nUnexpected under/overflow\n"); exit(1); }

	/* Calculate number of bits per character *************************/
	c=1; bits=0;
	do { c=c<<1; bits++; } while(c!=0);
	c= (char)(-1);
	printf("/\* Char = %d bits, %ssigned *\/\n", sizeof(c)*bits,
			((int)c)<0?"":"un");

	/* Shorts, ints and longs *****************************************/
	sprop();
	iprop();
	lprop();

	if (setjmp(lab)!=0) { printf("\nUnexpected under/overflow\n"); exit(1); }

	/* Alignment constants ********************************************/
	printf("/\* Alignments for char=%d short=%d int=%d long=%d *\/\n",
		sizeof(struct{char i; char c;})-sizeof(char),
		sizeof(struct{short i; char c;})-sizeof(short),
		sizeof(struct{int i; char c;})-sizeof(int),
		sizeof(struct{long i; char c;})-sizeof(long));

	/* Pointers *******************************************************/
	printf("/\* Char pointers = %d bits%s *\/\n", sizeof(char *)*bits,
		sizeof(char *)>sizeof(int)?" BEWARE! larger than int!":"");
	printf("/\* Int pointers = %d bits%s *\/\n", sizeof(int *)*bits,
		sizeof(int *)>sizeof(int)?" BEWARE! larger than int!":"");

	return bits;
}

int log(base, x) int base; double x; {
	int r=0;
	while (x>=base) { r++; x/=base; }
	return r;
}

int eprop() { /* See if expressions are evaluated in extended precision */
	int imant;
	double a, b, base;

	if (setjmp(lab)!=0) { printf("\nUnexpected under/overflow\n"); exit(1); }

	/* Size of mantissa **************************************/
	a=1.0;
	do { a=a+a; } while ((((a+1.0)-a)-1.0) == 0.0);
	b=1.0;
	do { b=b+b; } while ((base=((a+b)-a)) == 0.0);

	imant=0; b=1.0;
	do { imant++; b=b*base; }
	while ((((b+1.0)-b)-1.0) == 0.0);
	return imant;
}

#define fabs(x) (((x)<0.0)?(-x):(x))

#endif /* ifndef PASS */

/* As I said, I apologise for the contortions below. The procedures are
   expanded twice (for float and double) or three times (for short, int and
   long) by the preprocessor. That way, I never make a change to one that
   I forget to make to the other. #undef on an already undefined thing
   is (wrongly) flagged as an error by some compilers, therefore the #ifdef
   that follows: 
*/

#ifdef Number
#undef Number
#undef THING
#undef FPROP
#undef Store
#undef Sum
#undef Diff
#undef Mul
#undef Div
#endif

#ifdef Integer
#undef Integer
#undef INT
#undef IPROP
#endif

#if PASS == 1

#define Number float
#define THING "float"
#define FPROP fprop
#define Store fStore
#define Sum fSum
#define Diff fDiff
#define Mul fMul
#define Div fDiv

#define Integer short
#define INT "short"
#define IPROP sprop

#endif /* PASS == 1 */

#if PASS == 2

#define Number double
#define THING "double"
#define FPROP dprop
#define Store dStore
#define Sum dSum
#define Diff dDiff
#define Mul dMul
#define Div dDiv

#define Integer int
#define INT "int"
#define IPROP iprop

#endif /* if PASS == 2 */

#if PASS == 3

#define Integer long
#define INT "long"
#define IPROP lprop

#endif /* if PASS == 3 */

IPROP() {
	Integer newi, maxi, maxeri;
	int ibits, ipower, two=2;

	/* Calculate max short/int/long ***********************************/
	/* Calculate 2**n-1 until overflow - then use the previous value  */

	newi=1; maxi=0;

	if (setjmp(lab)==0)
		for(ipower=0; newi>maxi; ipower++) {
			maxi=newi;
			newi=newi*two+1;
		}

	/* Now for those daft Cybers: */

	maxeri=0; newi=maxi;

	if (setjmp(lab)==0)
		for(ibits=ipower; newi>maxeri; ibits++) {
			maxeri=newi;
			newi=newi+newi+1;
		}

	printf("/\* Maximum %s = %ld (= 2**%d-1) *\/\n",
				INT, (long)maxi, ipower);

	if (maxeri>maxi) {
		printf("/\* There is a larger %s, %ld (= 2**%d-1), %s *\/\n",
			INT, (long)maxeri, ibits, 
			"but only for addition, not multiplication");
	}
}

#ifdef Number

/* These routines are intended to defeat any attempt at optimisation */
Store(a, b) Number a, *b; { *b=a; }
Number Sum(a, b) Number a, b; { Number r; Store(a+b, &r); return (r); }
Number Diff(a, b) Number a, b; { Number r; Store(a-b, &r); return (r); }
Number Mul(a, b) Number a, b; { Number r; Store(a*b, &r); return (r); }
Number Div(a, b) Number a, b; { Number r; Store(a/b, &r); return (r); }

int FPROP(bits) int bits; {
	/* Properties of floating types, using algorithms by Cody and Waite
	   from MA Malcolm, as modified by WM Gentleman and SB Marovich.
	   Further extended by S Pemberton.

	   Returns the number of digits in the fraction.
	*/

	int i, ibase, iexp, irnd, imant, iz, k, machep, maxexp, minexp,
	    mx, negeps, mantbits;
	Number a, b, base, basein, basem1, eps, epsneg, xmax, newxmax,
	       xmin, xminner, y, y1, z, z1, z2;

	if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }

	printf("\n/\* Properties of %s: *\/\n", THING);

	/* Base and size of mantissa **************************************/
	a=1.0;
	do { a=Sum(a, a); } while (Diff(Diff(Sum(a, 1.0), a), 1.0) == 0.0);
	b=1.0;
	do { b=Sum(b, b); } while ((base=Diff(Sum(a, b), a)) == 0.0);
	ibase=base;
	printf("/\* Base = %d *\/\n", ibase);

	imant=0; b=1.0;
	do { imant++; b=Mul(b, base); }
	while (Diff(Diff(Sum(b,1.0),b),1.0) == 0.0);
	printf("/\* Significant base digits = %d %s %d %s *\/\n",
			imant, "(= at least", log(10, (double)b),
			"decimal digits)");

	/* Various flavours of epsilon ************************************/
	basem1=Diff(base,1.0);
	if (Diff(Sum(a, basem1), a) != 0.0) irnd=1; 
	else irnd=0;

	negeps=imant+imant;
	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;
	printf("/\* Smallest x such that 1.0-base**x != 1.0 = %d *\/\n", negeps);

	epsneg=a;
	if ((ibase!=2) && (irnd==1)) {
	/*	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;
	}
	printf("/\* Small x such that 1.0-x != 1.0 = %g *\/\n", epsneg);
	/* it may not be the smallest */

	machep= -imant-imant;
	a=b;
	while (Diff(Sum(1.0, a), 1.0) == 0.0) { a*=base; machep++; }
	printf("/\* Smallest x such that 1.0+base**x != 1.0 = %d *\/\n", machep);

	eps=a;
	if ((ibase!=2) && (irnd==1)) {
	/*	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) eps=a; => */
		if (Diff(Sum(1.0, a), 1.0) != 0.0) eps=a;
	}
	printf("/\* Smallest x such that 1.0+x != 1.0 = %g *\/\n", eps);

	/* Round or chop **************************************************/
	if (irnd == 1) { printf("/\* Arithmetic rounds *\/\n"); }
	else { 
		printf("/\* Arithmetic chops");
		if (Diff(Mul(Sum(1.0,eps),1.0),1.0) !=  0.0) {
			printf(" but uses guard digits");
		}
		printf(" *\/\n");
	}

	/* Size of and minimum normalised exponent ************************/
	y=0; i=0; k=1; z=basein; z1=(1.0+eps)/base;

	/* Coarse search for the largest power of two */
	if (setjmp(lab)==0) /* in case of underflow trap */
		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);

	if (ibase != 10) {
		iexp=i+1; /* for the sign */
		mx=k+k;
	} else {
		iexp=2;
		iz=ibase;
		while (k >= iz) { iz*=ibase; iexp++; }
		mx=iz+iz-1;
	}

	/* Fine tune starting with y and y1 */
	if (setjmp(lab)==0) /* in case of underflow trap */
		do {
			xmin=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) >= xmin)) break;
			k++;
		} while (1);

	if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }

	minexp=(-k)+1;

	if ((mx <= k+k-3) && (ibase != 10)) { mx+=mx; iexp+=1; }
	printf("/\* Number of bits used for exponent = %d *\/\n", iexp);
	printf("/\* Minimum normalised exponent = %d *\/\n", minexp);
	printf("/\* Minimum normalised positive number = %g *\/\n", xmin);

	/* Minimum exponent ************************************************/
	if (setjmp(lab)==0) /* in case of underflow trap */
		do {
			xminner=y;
			y=Div(y,base);
			a=Mul(y,1.0);
			if ((Sum(a,a) == 0.0) || (fabs(y) >= xminner)) break;
		} while (1);

	if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }

	if (xminner != 0.0 && xminner != xmin) {
		printf("/\* The smallest numbers are not kept normalised *\/\n");
		printf("/\* Smallest unnormalised positive number = %g *\/\n",
			xminner);
	} else printf("/\* The smallest numbers are normalised *\/\n");

#ifndef UNITY
	/* Maximum exponent ************************************************/
	maxexp=2; xmax=1.0; newxmax=base+1.0;
	if (setjmp(lab) == 0) {
		while (xmax<newxmax) {
			xmax=newxmax;
			newxmax=Mul(newxmax, base);
			if (Div(newxmax, base) != xmax) break; /* ieee infinity */
			maxexp++;
		}
	}
	if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }

	printf("/\* Maximum exponent = %d *\/\n", maxexp);

	/* Largest number ***************************************************/
	xmax=Diff(1.0, epsneg);
	if (Mul(xmax,1.0) != xmax) xmax=Diff(1.0, Mul(base,epsneg));
	for (i=1; i<=maxexp; i++) xmax=Mul(xmax, base);
	printf("/\* Maximum number = %g *\/\n", xmax);
#endif

	/* Hidden bit + sanity check ****************************************/
	if (ibase != 10) {
		mantbits=log(2, (double)ibase)*imant;
		if (mantbits+iexp+1 == sizeof(Number)*bits+1) {
			printf("/\* Arithmetic uses a hidden bit *\/\n");
		} else if (mantbits+iexp+1 == sizeof(Number)*bits) {
			printf("/\* Arithmetic doesn't use a hidden bit *\/\n");
		} else {
			printf("/\* Something fishy here! %s %s %s *\/\n",
				"Exponent size + mantissa size doesn't match",
				"with the size of a", THING);
		}
	}
	return imant;
}

#endif /* ifdef Number */

#if PASS == 3
#undef PASS
#define PASS 4
#endif

#if PASS == 2
#undef PASS
#define PASS 3
#endif

#if PASS == 1
#undef PASS
#define PASS 2
#endif

#if PASS < 4
#include __FILE__
#endif
DataMuseum.dk

DKUUG/EUUG Conference tapes

⟦cedff05a8⟧ TextFile

Derivation

TextFile
