/* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */ /* \$Header: mkconfig.c,v 1.4 85/08/26 10:41:39 timo Exp \$ */ /* Generate constants for configuration file */ /* If your C system is not unix but does have signal/setjmp, */ /* add a #define unix */ /* You may also need to add some calls to signal(). */ #ifdef unix #define SIGNAL #include #include jmp_buf lab; overflow(sig) int sig; { /* what to do on overflow/underflow */ signal(sig, overflow); longjmp(lab, 1); } #else /* Dummy routines instead */ int lab=1; int setjmp(lab) int lab; { return(0); } #endif #define fabs(x) (((x)<0.0)?(-x):(x)) #define min(x,y) (((x)<(y))?(x):(y)) /* These routines are intended to defeat any attempt at optimisation */ Dstore(a, b) double a, *b; { *b=a; } double Dsum(a, b) double a, b; { double r; Dstore(a+b, &r); return (r); } double Ddiff(a, b) double a, b; { double r; Dstore(a-b, &r); return (r); } double Dmul(a, b) double a, b; { double r; Dstore(a*b, &r); return (r); } double Ddiv(a, b) double a, b; { double r; Dstore(a/b, &r); return (r); } double power(x, n) int x, n; { double r=1.0; for (;n>0; n--) r*=x; return r; } int log(base, x) int base; double x; { int r=0; while (x>=base) { r++; x/=base; } return r; } main(argc, argv) int argc; char *argv[]; { char c; short newshort, maxshort, maxershort; int newint, maxint, shortbits, bits, mantbits, *p, shortpower, intpower, longpower; long newlong, maxlong, count; int i, ibase, iexp, irnd, imant, iz, k, machep, maxexp, minexp, mx, negeps, ngrd, normalised; double a, b, base, basein, basem1, eps, epsneg, xmax, newxmax, xmin, xminner, y, y1, z, z1, z2; double BIG, Maxreal; int BASE, MAXNUMDIG, ipower, tenlogBASE, Maxexpo, Minexpo, Dblbits; #ifdef SIGNAL signal(SIGFPE, overflow); if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); } #endif /****** Calculate max short *********************************************/ /* Calculate 2**n-1 until overflow - then use the previous value */ newshort=1; maxshort=0; if (setjmp(lab)==0) for(shortpower=0; newshort>maxshort; shortpower++) { maxshort=newshort; newshort=newshort*2+1; } /* Now for those daft Cybers: */ maxershort=0; newshort=maxshort; if (setjmp(lab)==0) for(shortbits=shortpower; newshort>maxershort; shortbits++) { maxershort=newshort; newshort=newshort+newshort+1; } bits= (shortbits+1)/sizeof(short); c= (char)(-1); printf("/\* char=%d bits, %ssigned *\/\n", sizeof(c)*bits, ((int)c)<0?"":"un"); printf("/\* maxshort=%d (=2**%d-1) *\/\n", maxshort, shortpower); if (maxershort>maxshort) { printf("/\* There is a larger maxshort, %d (=2**%d-1), %s *\/\n", maxershort, shortbits, "but only for addition, not multiplication"); } /****** Calculate max int by the same method ***************************/ newint=1; maxint=0; if (setjmp(lab)==0) for(intpower=0; newint>maxint; intpower++) { maxint=newint; newint=newint*2+1; } printf("/\* maxint=%d (=2**%d-1) *\/\n", maxint, intpower); /****** Calculate max long by the same method ***************************/ newlong=1; maxlong=0; if (setjmp(lab)==0) for(longpower=0; newlong>maxlong; longpower++) { maxlong=newlong; newlong=newlong*2+1; } if (setjmp(lab)!=0) { printf("\nUnexpected under/overflow\n"); exit(1); } printf("/\* maxlong=%ld (=2**%d-1) *\/\n", maxlong, longpower); /****** Pointers ********************************************************/ printf("/\* pointers=%d bits%s *\/\n", sizeof(p)*bits, sizeof(p)>sizeof(int)?" BEWARE! larger than int!":""); /****** Base and size of mantissa ***************************************/ a=1.0; do { a=Dsum(a, a); } while (Ddiff(Ddiff(Dsum(a, 1.0), a), 1.0) == 0.0); b=1.0; do { b=Dsum(b, b); } while ((base=Ddiff(Dsum(a, b), a)) == 0.0); ibase=base; printf("/\* base=%d *\/\n", ibase); imant=0; b=1.0; do { imant++; b=Dmul(b, base); } while (Ddiff(Ddiff(Dsum(b,1.0),b),1.0) == 0.0); printf("/\* Significant base digits=%d *\/\n", imant); /****** Various flavours of epsilon *************************************/ basem1=Ddiff(base,1.0); if (Ddiff(Dsum(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 (Ddiff(Ddiff(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=Ddiv(Dmul(a, Dsum(1.0, a)), Dsum(1.0, 1.0)); /* if ((1.0-a)-1.0 != 0.0) epsneg=a; => */ if (Ddiff(Ddiff(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 (Ddiff(Dsum(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=Ddiv(Dmul(a, Dsum(1.0, a)), Dsum(1.0, 1.0)); /* if ((1.0+a)-1.0 != 0.0) eps=a; => */ if (Ddiff(Dsum(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 ***************************************************/ ngrd=0; if (irnd == 1) { printf("/\* Arithmetic rounds *\/\n"); } else { printf("/\* Arithmetic chops"); if (Ddiff(Dmul(Dsum(1.0,eps),1.0),1.0) != 0.0) { ngrd=1; 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=Dmul(y,y); z1=Dmul(z1, y); a=Dmul(z,1.0); z2=Ddiv(z1,y); if (z2 != y1) break; if ((Dsum(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=Ddiv(y,base); y1=Ddiv(y1,base); a=Dmul(y,1.0); z2=Dmul(y1,base); if (z2 != z1) break; if ((Dsum(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=Ddiv(y,base); a=Dmul(y,1.0); if ((Dsum(a,a) == 0.0) || (fabs(y) >= xminner)) break; } while (1); if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); } normalised=1; if (xminner != 0.0 && xminner != xmin) { normalised=0; printf("/\* The smallest numbers are not kept normalised *\/\n"); printf("/\* Smallest unnormalised positive number=%g *\/\n", xminner); } /****** Maximum exponent ***************************************************/ maxexp=2; xmax=1.0; newxmax=base+1.0; if (setjmp(lab) == 0) { while (xmax