4.3BSD - /usr/src/new/B/src/bint/mkconfig.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: mkconfig.c,v 1.4 85/08/26 10:41:39 timo Exp $
   5: */
   6: 
   7: /* Generate constants for configuration file			*/
   8: 
   9: /* If your C system is not unix but does have signal/setjmp,	*/
  10: /*    add a #define unix					*/
  11: /* You may also need to add some calls to signal().		*/
  12: 
  13: #ifdef unix
  14: 
  15: #define SIGNAL
  16: 
  17: #include <signal.h>
  18: #include <setjmp.h>
  19: 
  20:     jmp_buf lab;
  21:     overflow(sig) int sig; { /* what to do on overflow/underflow */
  22:         signal(sig, overflow);
  23:         longjmp(lab, 1);
  24:     }
  25: 
  26: #else
  27:     /* Dummy routines instead */
  28:     int lab=1;
  29:     int setjmp(lab) int lab; { return(0); }
  30: 
  31: #endif
  32: 
  33: #define fabs(x) (((x)<0.0)?(-x):(x))
  34: #define min(x,y) (((x)<(y))?(x):(y))
  35: 
  36: /* These routines are intended to defeat any attempt at optimisation */
  37: Dstore(a, b) double a, *b; { *b=a; }
  38: double Dsum(a, b) double a, b; { double r; Dstore(a+b, &r); return (r); }
  39: double Ddiff(a, b) double a, b; { double r; Dstore(a-b, &r); return (r); }
  40: double Dmul(a, b) double a, b; { double r; Dstore(a*b, &r); return (r); }
  41: double Ddiv(a, b) double a, b; { double r; Dstore(a/b, &r); return (r); }
  42: 
  43: double power(x, n) int x, n; {
  44:     double r=1.0;
  45:     for (;n>0; n--) r*=x;
  46:     return r;
  47: }
  48: 
  49: int log(base, x) int base; double x; {
  50:     int r=0;
  51:     while (x>=base) { r++; x/=base; }
  52:     return r;
  53: }
  54: 
  55: main(argc, argv) int argc; char *argv[]; {
  56:     char c;
  57:     short newshort, maxshort, maxershort;
  58:     int newint, maxint, shortbits, bits, mantbits,
  59:         *p, shortpower, intpower, longpower;
  60:     long newlong, maxlong, count;
  61: 
  62:     int i, ibase, iexp, irnd, imant, iz, k, machep, maxexp, minexp,
  63:         mx, negeps, ngrd, normalised;
  64:     double a, b, base, basein, basem1, eps, epsneg, xmax, newxmax,
  65:            xmin, xminner, y, y1, z, z1, z2;
  66: 
  67:     double BIG, Maxreal;
  68:     int BASE, MAXNUMDIG, ipower, tenlogBASE, Maxexpo, Minexpo, Dblbits;
  69: 
  70: #ifdef SIGNAL
  71:     signal(SIGFPE, overflow);
  72:     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
  73: #endif
  74: 
  75: /****** Calculate max short *********************************************/
  76: /*      Calculate 2**n-1 until overflow - then use the previous value	*/
  77: 
  78:     newshort=1; maxshort=0;
  79: 
  80:     if (setjmp(lab)==0)
  81:         for(shortpower=0; newshort>maxshort; shortpower++) {
  82:             maxshort=newshort;
  83:             newshort=newshort*2+1;
  84:         }
  85: 
  86:     /* Now for those daft Cybers: */
  87: 
  88:     maxershort=0; newshort=maxshort;
  89: 
  90:     if (setjmp(lab)==0)
  91:         for(shortbits=shortpower; newshort>maxershort; shortbits++) {
  92:             maxershort=newshort;
  93:             newshort=newshort+newshort+1;
  94:         }
  95: 
  96:     bits= (shortbits+1)/sizeof(short);
  97:     c= (char)(-1);
  98:     printf("/\* char=%d bits, %ssigned *\/\n", sizeof(c)*bits,
  99:             ((int)c)<0?"":"un");
 100:     printf("/\* maxshort=%d (=2**%d-1) *\/\n", maxshort, shortpower);
 101: 
 102:     if (maxershort>maxshort) {
 103:         printf("/\* There is a larger maxshort, %d (=2**%d-1), %s *\/\n",
 104:             maxershort, shortbits,
 105:             "but only for addition, not multiplication");
 106:     }
 107: 
 108: /****** Calculate max int by the same method ***************************/
 109: 
 110:     newint=1; maxint=0;
 111: 
 112:     if (setjmp(lab)==0)
 113:         for(intpower=0; newint>maxint; intpower++) {
 114:             maxint=newint;
 115:             newint=newint*2+1;
 116:         }
 117: 
 118:     printf("/\* maxint=%d (=2**%d-1) *\/\n", maxint, intpower);
 119: 
 120: /****** Calculate max long by the same method ***************************/
 121: 
 122:     newlong=1; maxlong=0;
 123: 
 124:     if (setjmp(lab)==0)
 125:         for(longpower=0; newlong>maxlong; longpower++) {
 126:             maxlong=newlong;
 127:             newlong=newlong*2+1;
 128:         }
 129: 
 130:     if (setjmp(lab)!=0) { printf("\nUnexpected under/overflow\n"); exit(1); }
 131: 
 132:     printf("/\* maxlong=%ld (=2**%d-1) *\/\n", maxlong, longpower);
 133: 
 134: /****** Pointers ********************************************************/
 135:     printf("/\* pointers=%d bits%s *\/\n", sizeof(p)*bits,
 136:         sizeof(p)>sizeof(int)?" BEWARE! larger than int!":"");
 137: 
 138: /****** Base and size of mantissa ***************************************/
 139:     a=1.0;
 140:     do { a=Dsum(a, a); } while (Ddiff(Ddiff(Dsum(a, 1.0), a), 1.0) == 0.0);
 141:     b=1.0;
 142:     do { b=Dsum(b, b); } while ((base=Ddiff(Dsum(a, b), a)) == 0.0);
 143:     ibase=base;
 144:     printf("/\* base=%d *\/\n", ibase);
 145: 
 146:     imant=0; b=1.0;
 147:     do { imant++; b=Dmul(b, base); }
 148:     while (Ddiff(Ddiff(Dsum(b,1.0),b),1.0) == 0.0);
 149:     printf("/\* Significant base digits=%d *\/\n", imant);
 150: 
 151: /****** Various flavours of epsilon *************************************/
 152:     basem1=Ddiff(base,1.0);
 153:     if (Ddiff(Dsum(a, basem1), a) != 0.0) irnd=1;
 154:     else irnd=0;
 155: 
 156:     negeps=imant+imant;
 157:     basein=1.0/base;
 158:     a=1.0;
 159:     for(i=1; i<=negeps; i++) a*=basein;
 160: 
 161:     b=a;
 162:     while (Ddiff(Ddiff(1.0, a), 1.0) == 0.0) {
 163:         a*=base;
 164:         negeps--;
 165:     }
 166:     negeps= -negeps;
 167:     printf("/\* Smallest x such that 1.0-base**x != 1.0=%d *\/\n", negeps);
 168: 
 169:     epsneg=a;
 170:     if ((ibase!=2) && (irnd==1)) {
 171:     /*	a=(a*(1.0+a))/(1.0+1.0); => */
 172:         a=Ddiv(Dmul(a, Dsum(1.0, a)), Dsum(1.0, 1.0));
 173:     /*	if ((1.0-a)-1.0 != 0.0) epsneg=a; => */
 174:         if (Ddiff(Ddiff(1.0, a), 1.0) != 0.0) epsneg=a;
 175:     }
 176:     printf("/\* Small x such that 1.0-x != 1.0=%g *\/\n", epsneg);
 177:     /* it may not be the smallest */
 178: 
 179:     machep= -imant-imant;
 180:     a=b;
 181:     while (Ddiff(Dsum(1.0, a), 1.0) == 0.0) { a*=base; machep++; }
 182:     printf("/\* Smallest x such that 1.0+base**x != 1.0=%d *\/\n", machep);
 183: 
 184:     eps=a;
 185:     if ((ibase!=2) && (irnd==1)) {
 186:     /*	a=(a*(1.0+a))/(1.0+1.0); => */
 187:         a=Ddiv(Dmul(a, Dsum(1.0, a)), Dsum(1.0, 1.0));
 188:     /*	if ((1.0+a)-1.0 != 0.0) eps=a; => */
 189:         if (Ddiff(Dsum(1.0, a), 1.0) != 0.0) eps=a;
 190:     }
 191:     printf("/\* Smallest x such that 1.0+x != 1.0=%g *\/\n", eps);
 192: 
 193: /****** Round or chop ***************************************************/
 194:     ngrd=0;
 195:     if (irnd == 1) { printf("/\* Arithmetic rounds *\/\n"); }
 196:     else {
 197:         printf("/\* Arithmetic chops");
 198:         if (Ddiff(Dmul(Dsum(1.0,eps),1.0),1.0) != 0.0) {
 199:             ngrd=1;
 200:             printf(" but uses guard digits");
 201:         }
 202:         printf(" *\/\n");
 203:     }
 204: 
 205: /****** Size of and minimum normalised exponent ****************************/
 206:     y=0; i=0; k=1; z=basein; z1=(1.0+eps)/base;
 207: 
 208:     /* Coarse search for the largest power of two */
 209:     if (setjmp(lab)==0) /* in case of underflow trap */
 210:         do {
 211:             y=z; y1=z1;
 212:             z=Dmul(y,y); z1=Dmul(z1, y);
 213:             a=Dmul(z,1.0);
 214:             z2=Ddiv(z1,y);
 215:             if (z2 != y1) break;
 216:             if ((Dsum(a,a) == 0.0) || (fabs(z) >= y)) break;
 217:             i++;
 218:             k+=k;
 219:         } while(1);
 220: 
 221:     if (ibase != 10) {
 222:         iexp=i+1; /* for the sign */
 223:         mx=k+k;
 224:     } else {
 225:         iexp=2;
 226:         iz=ibase;
 227:         while (k >= iz) { iz*=ibase; iexp++; }
 228:         mx=iz+iz-1;
 229:     }
 230: 
 231:     /* Fine tune starting with y and y1 */
 232:     if (setjmp(lab)==0) /* in case of underflow trap */
 233:         do {
 234:             xmin=y; z1=y1;
 235:             y=Ddiv(y,base); y1=Ddiv(y1,base);
 236:             a=Dmul(y,1.0);
 237:             z2=Dmul(y1,base);
 238:             if (z2 != z1) break;
 239:             if ((Dsum(a,a) == 0.0) || (fabs(y) >= xmin)) break;
 240:             k++;
 241:         } while (1);
 242: 
 243:     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
 244: 
 245:     minexp= (-k)+1;
 246: 
 247:     if ((mx <= k+k-3) && (ibase != 10)) { mx+=mx; iexp+=1; }
 248:     printf("/\* Number of bits used for exponent=%d *\/\n", iexp);
 249:     printf("/\* Minimum normalised exponent=%d *\/\n", minexp);
 250:     printf("/\* Minimum normalised positive number=%g *\/\n", xmin);
 251: 
 252: /****** Minimum exponent ***************************************************/
 253:     if (setjmp(lab)==0) /* in case of underflow trap */
 254:         do {
 255:             xminner=y;
 256:             y=Ddiv(y,base);
 257:             a=Dmul(y,1.0);
 258:             if ((Dsum(a,a) == 0.0) || (fabs(y) >= xminner)) break;
 259:         } while (1);
 260: 
 261:     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
 262: 
 263:     normalised=1;
 264:     if (xminner != 0.0 && xminner != xmin) {
 265:         normalised=0;
 266:         printf("/\* The smallest numbers are not kept normalised *\/\n");
 267:         printf("/\* Smallest unnormalised positive number=%g *\/\n",
 268:             xminner);
 269:     }
 270: 
 271: /****** Maximum exponent ***************************************************/
 272:     maxexp=2; xmax=1.0; newxmax=base+1.0;
 273:     if (setjmp(lab) == 0) {
 274:         while (xmax<newxmax) {
 275:             xmax=newxmax;
 276:             newxmax=Dmul(newxmax, base);
 277:             if (Ddiv(newxmax, base) != xmax) break; /* ieee infinity */
 278:             maxexp++;
 279:         }
 280:     }
 281:     if (setjmp(lab)!=0) { printf("Unexpected over/underflow\n"); exit(1); }
 282: 
 283:     printf("/\* Maximum exponent=%d *\/\n", maxexp);
 284: 
 285: /****** Largest and smallest numbers ************************************/
 286:     xmax=Ddiff(1.0, epsneg);
 287:     if (Dmul(xmax,1.0) != xmax) xmax=Ddiff(1.0, Dmul(base,epsneg));
 288:     for (i=1; i<=maxexp; i++) xmax=Dmul(xmax, base);
 289:     printf("/\* Maximum number=%g *\/\n", xmax);
 290: 
 291: /****** Hidden bit + sanity check ***************************************/
 292:     if (ibase != 10) {
 293:         mantbits=log(2, (double)ibase)*imant;
 294:         if (mantbits+iexp+1 == sizeof(double)*bits+1) {
 295:             printf("/\* Double arithmetic uses a hidden bit *\/\n");
 296:         } else if (mantbits+iexp+1 == sizeof(double)*bits) {
 297:             printf("/\* Double arithmetic doesn't use a hidden bit *\/\n");
 298:         } else {
 299:             printf("/\* Something fishy here! %s %s *\/\n",
 300:                 "Exponent size + mantissa size doesn't match",
 301:                 "with the size of a double.");
 302:         }
 303:     }
 304: 
 305: /****** The point of it all: ********************************************/
 306:     printf("\n/\* Numeric package constants *\/\n");
 307: 
 308:     BIG= power(ibase, imant)-1.0;
 309:     MAXNUMDIG= log(10, BIG);
 310:     ipower= log(10, (double)(maxint/2));
 311:     Maxreal= xmax;
 312:     Maxexpo= log(2, (double)ibase)*maxexp;
 313:     Minexpo= log(2, (double)ibase)*minexp;
 314:     Dblbits= log(2, (double)ibase)*imant;
 315:     tenlogBASE= min(MAXNUMDIG/2, ipower);
 316:     BASE=1; for(i=1; i<=tenlogBASE; i++) BASE*=10;
 317: 
 318:     printf("#define BIG %1.1f /\* Maximum integral double *\/\n", BIG);
 319:     printf("#define LONG ");
 320:         for(i=1; i<=MAXNUMDIG; i++) printf("9");
 321:         printf(".5 /\* Maximum power of ten less than BIG *\/\n");
 322:     printf("#define MAXNUMDIG %d /\* The number of 9's in LONG *\/\n",
 323:         MAXNUMDIG);
 324:     printf("#define MINNUMDIG 6 /\* Don't change: this is here for consistency *\/\n");
 325: 
 326:     printf("#define BASE %d /\* %s *\/\n",
 327:         BASE, "BASE**2<=BIG && BASE*2<=Maxint && -2*BASE>=Minint");
 328:     if (((double)BASE)*BASE > BIG || ((double)BASE)+BASE > maxint) {
 329:         printf("*** BASE value wrong\n");
 330:         exit(1);
 331:     }
 332:     printf("#define tenlogBASE %d /\*  = log10(BASE) *\/\n", tenlogBASE);
 333:     printf("#define Maxreal %e /\* Maximum double *\/\n", Maxreal);
 334:     printf("#define Maxexpo %d /\* Maximum value such that 2**Maxexpo<=Maxreal *\/\n",
 335:         Maxexpo);
 336:     printf("#define Minexpo (%d) /\* Minimum value such that -2**Minexpo>=Minreal *\/\n",
 337:         Minexpo);
 338:     printf("#define Dblbits %d /\* Maximum value such that 2**Dblbits<=BIG *\/\n",
 339:         Dblbits);
 340: 
 341:     printf("#define Maxintlet %d /\* Maximum short *\/\n", maxshort);
 342:     printf("#define Maxint %d /\* Maximum int *\/\n", maxint);
 343: 
 344:     printf("#define Maxsmallint %d /\* BASE-1 *\/\n", BASE-1);
 345:     printf("#define Minsmallint (%d) /\* -BASE *\/\n", -BASE);
 346: 
 347: /* An extra goody: the approximate amount of data-space */
 348: /* Put here because it is likely to be slower then the rest */
 349: 
 350:     /*Allocate blocks of 1000 until no more available*/
 351:     /*Don't be tempted to change this to 1024: */
 352:     /*we don't know how much header information there is*/
 353: 
 354:     for(count=0; (p=(int *)malloc(1000))!=0; count++) { }
 355: 
 356:     printf("\n/\* Memory~= %d000 *\/\n", count);
 357: 
 358:     exit(0);
 359: }

Defined functions

Defined variables

Defined macros