```   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

Ddiff defined in line 39; used 16 times
Ddiv defined in line 41; used 7 times
Dmul defined in line 40; used 14 times
Dstore defined in line 37; used 4 times
Dsum defined in line 38; used 16 times
log defined in line 49; used 11 times
main defined in line 55; never used
overflow defined in line 21; used 2 times
power defined in line 43; used 1 times
setjmp defined in line 29; used 15 times

#### Defined variables

lab defined in line 28; used 16 times

#### Defined macros

SIGNAL defined in line 15; used 1 times
• in line 70
fabs defined in line 33; used 3 times
min defined in line 34; used 1 times
 Last modified: 1985-08-27 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 2666