![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to mkconfig.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint|
Return to mkconfig.c CVS log | Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint |
1.1 root 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.