|
709
|
1 /* |
|
|
2 |
|
|
3 This file combines the single and double precision versions of machar, |
|
|
4 selected by cc -DSP or cc -DDP. This feature provided by D. G. Hough, |
|
|
5 August 3, 1988. |
|
|
6 |
|
|
7 */ |
|
|
8 |
|
|
9 #ifdef SP |
|
|
10 #define REAL float |
|
|
11 #define ZERO 0.0 |
|
|
12 #define ONE 1.0 |
|
|
13 #define PREC "Single " |
|
|
14 #define REALSIZE 1 |
|
|
15 #endif |
|
|
16 |
|
|
17 #ifdef DP |
|
|
18 #define REAL double |
|
|
19 #define ZERO 0.0e0 |
|
|
20 #define ONE 1.0e0 |
|
|
21 #define PREC "Double " |
|
|
22 #define REALSIZE 2 |
|
|
23 #endif |
|
|
24 |
|
|
25 #include <math.h> |
|
|
26 #include <stdio.h> |
|
|
27 |
|
|
28 #define ABS(xxx) ((xxx>ZERO)?(xxx):(-xxx)) |
|
|
29 |
|
|
30 void |
|
|
31 rmachar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp, |
|
|
32 maxexp,eps,epsneg,xmin,xmax) |
|
|
33 |
|
|
34 int *ibeta,*iexp,*irnd,*it,*machep,*maxexp,*minexp,*negep,*ngrd; |
|
|
35 REAL *eps,*epsneg,*xmax,*xmin; |
|
|
36 |
|
|
37 /* |
|
|
38 |
|
|
39 This subroutine is intended to determine the parameters of the |
|
|
40 floating-point arithmetic system specified below. The |
|
|
41 determination of the first three uses an extension of an algorithm |
|
|
42 due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some, |
|
|
43 but not all, of the improvements suggested by M. Gentleman and S. |
|
|
44 Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this |
|
|
45 program was published in the book Software Manual for the |
|
|
46 Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall, |
|
|
47 Englewood Cliffs, NJ, 1980. The present program is a |
|
|
48 translation of the Fortran 77 program in W. J. Cody, "MACHAR: |
|
|
49 A subroutine to dynamically determine machine parameters". |
|
|
50 TOMS (14), 1988. |
|
|
51 |
|
|
52 Parameter values reported are as follows: |
|
|
53 |
|
|
54 ibeta - the radix for the floating-point representation |
|
|
55 it - the number of base ibeta digits in the floating-point |
|
|
56 significand |
|
|
57 irnd - 0 if floating-point addition chops |
|
|
58 1 if floating-point addition rounds, but not in the |
|
|
59 IEEE style |
|
|
60 2 if floating-point addition rounds in the IEEE style |
|
|
61 3 if floating-point addition chops, and there is |
|
|
62 partial underflow |
|
|
63 4 if floating-point addition rounds, but not in the |
|
|
64 IEEE style, and there is partial underflow |
|
|
65 5 if floating-point addition rounds in the IEEE style, |
|
|
66 and there is partial underflow |
|
|
67 ngrd - the number of guard digits for multiplication with |
|
|
68 truncating arithmetic. It is |
|
|
69 0 if floating-point arithmetic rounds, or if it |
|
|
70 truncates and only it base ibeta digits |
|
|
71 participate in the post-normalization shift of the |
|
|
72 floating-point significand in multiplication; |
|
|
73 1 if floating-point arithmetic truncates and more |
|
|
74 than it base ibeta digits participate in the |
|
|
75 post-normalization shift of the floating-point |
|
|
76 significand in multiplication. |
|
|
77 machep - the largest negative integer such that |
|
|
78 1.0+FLOAT(ibeta)**machep .NE. 1.0, except that |
|
|
79 machep is bounded below by -(it+3) |
|
|
80 negeps - the largest negative integer such that |
|
|
81 1.0-FLOAT(ibeta)**negeps .NE. 1.0, except that |
|
|
82 negeps is bounded below by -(it+3) |
|
|
83 iexp - the number of bits (decimal places if ibeta = 10) |
|
|
84 reserved for the representation of the exponent |
|
|
85 (including the bias or sign) of a floating-point |
|
|
86 number |
|
|
87 minexp - the largest in magnitude negative integer such that |
|
|
88 FLOAT(ibeta)**minexp is positive and normalized |
|
|
89 maxexp - the smallest positive power of BETA that overflows |
|
|
90 eps - the smallest positive floating-point number such |
|
|
91 that 1.0+eps .NE. 1.0. In particular, if either |
|
|
92 ibeta = 2 or IRND = 0, eps = FLOAT(ibeta)**machep. |
|
|
93 Otherwise, eps = (FLOAT(ibeta)**machep)/2 |
|
|
94 epsneg - A small positive floating-point number such that |
|
|
95 1.0-epsneg .NE. 1.0. In particular, if ibeta = 2 |
|
|
96 or IRND = 0, epsneg = FLOAT(ibeta)**negeps. |
|
|
97 Otherwise, epsneg = (ibeta**negeps)/2. Because |
|
|
98 negeps is bounded below by -(it+3), epsneg may not |
|
|
99 be the smallest number that can alter 1.0 by |
|
|
100 subtraction. |
|
|
101 xmin - the smallest non-vanishing normalized floating-point |
|
|
102 power of the radix, i.e., xmin = FLOAT(ibeta)**minexp |
|
|
103 xmax - the largest finite floating-point number. In |
|
|
104 particular xmax = (1.0-epsneg)*FLOAT(ibeta)**maxexp |
|
|
105 Note - on some machines xmax will be only the |
|
|
106 second, or perhaps third, largest number, being |
|
|
107 too small by 1 or 2 units in the last digit of |
|
|
108 the significand. |
|
|
109 |
|
|
110 Latest revision - August 4, 1988 |
|
|
111 |
|
|
112 Author - W. J. Cody |
|
|
113 Argonne National Laboratory |
|
|
114 |
|
|
115 */ |
|
|
116 |
|
|
117 { |
|
|
118 int i,iz,j,k; |
|
|
119 int mx,itmp,nxres; |
|
|
120 REAL a,b,beta,betain,one,y,z,zero; |
|
|
121 REAL betah,t,tmp,tmpa,tmp1,two; |
|
|
122 |
|
|
123 (*irnd) = 1; |
|
|
124 one = (REAL)(*irnd); |
|
|
125 two = one + one; |
|
|
126 a = two; |
|
|
127 b = a; |
|
|
128 zero = 0.0e0; |
|
|
129 |
|
|
130 /* |
|
|
131 determine ibeta,beta ala malcolm |
|
|
132 */ |
|
|
133 |
|
|
134 tmp = ((a+one)-a)-one; |
|
|
135 |
|
|
136 while (tmp == zero) { |
|
|
137 a = a+a; |
|
|
138 tmp = a+one; |
|
|
139 tmp1 = tmp-a; |
|
|
140 tmp = tmp1-one; |
|
|
141 } |
|
|
142 |
|
|
143 tmp = a+b; |
|
|
144 itmp = (int)(tmp-a); |
|
|
145 while (itmp == 0) { |
|
|
146 b = b+b; |
|
|
147 tmp = a+b; |
|
|
148 itmp = (int)(tmp-a); |
|
|
149 } |
|
|
150 |
|
|
151 *ibeta = itmp; |
|
|
152 beta = (REAL)(*ibeta); |
|
|
153 |
|
|
154 /* |
|
|
155 determine irnd, it |
|
|
156 */ |
|
|
157 |
|
|
158 (*it) = 0; |
|
|
159 b = one; |
|
|
160 tmp = ((b+one)-b)-one; |
|
|
161 |
|
|
162 while (tmp == zero) { |
|
|
163 *it = *it+1; |
|
|
164 b = b*beta; |
|
|
165 tmp = b+one; |
|
|
166 tmp1 = tmp-b; |
|
|
167 tmp = tmp1-one; |
|
|
168 } |
|
|
169 |
|
|
170 *irnd = 0; |
|
|
171 betah = beta/two; |
|
|
172 tmp = a+betah; |
|
|
173 tmp1 = tmp-a; |
|
|
174 if (tmp1 != zero) *irnd = 1; |
|
|
175 tmpa = a+beta; |
|
|
176 tmp = tmpa+betah; |
|
|
177 if ((*irnd == 0) && (tmp-tmpa != zero)) *irnd = 2; |
|
|
178 |
|
|
179 /* |
|
|
180 determine negep, epsneg |
|
|
181 */ |
|
|
182 |
|
|
183 (*negep) = (*it) + 3; |
|
|
184 betain = one / beta; |
|
|
185 a = one; |
|
|
186 |
|
|
187 for (i = 1; i<=(*negep); i++) { |
|
|
188 a = a * betain; |
|
|
189 } |
|
|
190 |
|
|
191 b = a; |
|
|
192 tmp = (one-a); |
|
|
193 tmp = tmp-one; |
|
|
194 |
|
|
195 while (tmp == zero) { |
|
|
196 a = a*beta; |
|
|
197 *negep = *negep-1; |
|
|
198 tmp1 = one-a; |
|
|
199 tmp = tmp1-one; |
|
|
200 } |
|
|
201 |
|
|
202 (*negep) = -(*negep); |
|
|
203 (*epsneg) = a; |
|
|
204 |
|
|
205 /* |
|
|
206 determine machep, eps |
|
|
207 */ |
|
|
208 |
|
|
209 (*machep) = -(*it) - 3; |
|
|
210 a = b; |
|
|
211 tmp = one+a; |
|
|
212 |
|
|
213 while (tmp-one == zero) { |
|
|
214 a = a*beta; |
|
|
215 *machep = *machep+1; |
|
|
216 tmp = one+a; |
|
|
217 } |
|
|
218 |
|
|
219 *eps = a; |
|
|
220 |
|
|
221 /* |
|
|
222 determine ngrd |
|
|
223 */ |
|
|
224 |
|
|
225 (*ngrd) = 0; |
|
|
226 tmp = one+*eps; |
|
|
227 tmp = tmp*one; |
|
|
228 if (((*irnd) == 0) && (tmp-one) != zero) (*ngrd) = 1; |
|
|
229 |
|
|
230 /* |
|
|
231 determine iexp, minexp, xmin |
|
|
232 |
|
|
233 loop to determine largest i such that |
|
|
234 (1/beta) ** (2**(i)) |
|
|
235 does not underflow. |
|
|
236 exit from loop is signaled by an underflow. |
|
|
237 */ |
|
|
238 |
|
|
239 i = 0; |
|
|
240 k = 1; |
|
|
241 z = betain; |
|
|
242 t = one+*eps; |
|
|
243 nxres = 0; |
|
|
244 |
|
|
245 for (;;) { |
|
|
246 y = z; |
|
|
247 z = y * y; |
|
|
248 |
|
|
249 /* |
|
|
250 check for underflow |
|
|
251 */ |
|
|
252 |
|
|
253 a = z * one; |
|
|
254 tmp = z*t; |
|
|
255 if ((a+a == zero) || (ABS(z) > y)) break; |
|
|
256 tmp1 = tmp*betain; |
|
|
257 if (tmp1*beta == z) break; |
|
|
258 i = i + 1; |
|
|
259 k = k+k; |
|
|
260 } |
|
|
261 |
|
|
262 /* |
|
|
263 determine k such that (1/beta)**k does not underflow |
|
|
264 first set k = 2 ** i |
|
|
265 */ |
|
|
266 |
|
|
267 (*iexp) = i + 1; |
|
|
268 mx = k + k; |
|
|
269 if (*ibeta == 10) { |
|
|
270 |
|
|
271 /* |
|
|
272 for decimal machines only |
|
|
273 */ |
|
|
274 |
|
|
275 (*iexp) = 2; |
|
|
276 iz = *ibeta; |
|
|
277 while (k >= iz) { |
|
|
278 iz = iz * (*ibeta); |
|
|
279 (*iexp) = (*iexp) + 1; |
|
|
280 } |
|
|
281 mx = iz + iz - 1; |
|
|
282 } |
|
|
283 |
|
|
284 /* |
|
|
285 loop to determine minexp, xmin. |
|
|
286 exit from loop is signaled by an underflow. |
|
|
287 */ |
|
|
288 |
|
|
289 for (;;) { |
|
|
290 (*xmin) = y; |
|
|
291 y = y * betain; |
|
|
292 a = y * one; |
|
|
293 tmp = y*t; |
|
|
294 tmp1 = a+a; |
|
|
295 if ((tmp1 == zero) || (ABS(y) >= (*xmin))) break; |
|
|
296 k = k + 1; |
|
|
297 tmp1 = tmp*betain; |
|
|
298 tmp1 = tmp1*beta; |
|
|
299 |
|
|
300 if ((tmp1 == y) && (tmp != y)) { |
|
|
301 nxres = 3; |
|
|
302 *xmin = y; |
|
|
303 break; |
|
|
304 } |
|
|
305 |
|
|
306 } |
|
|
307 |
|
|
308 (*minexp) = -k; |
|
|
309 |
|
|
310 /* |
|
|
311 determine maxexp, xmax |
|
|
312 */ |
|
|
313 |
|
|
314 if ((mx <= k+k-3) && ((*ibeta) != 10)) { |
|
|
315 mx = mx + mx; |
|
|
316 (*iexp) = (*iexp) + 1; |
|
|
317 } |
|
|
318 |
|
|
319 (*maxexp) = mx + (*minexp); |
|
|
320 |
|
|
321 /* |
|
|
322 Adjust *irnd to reflect partial underflow. |
|
|
323 */ |
|
|
324 |
|
|
325 (*irnd) = (*irnd)+nxres; |
|
|
326 |
|
|
327 /* |
|
|
328 Adjust for IEEE style machines. |
|
|
329 */ |
|
|
330 |
|
|
331 if ((*irnd) >= 2) (*maxexp) = (*maxexp)-2; |
|
|
332 |
|
|
333 /* |
|
|
334 adjust for machines with implicit leading bit in binary |
|
|
335 significand and machines with radix point at extreme |
|
|
336 right of significand. |
|
|
337 */ |
|
|
338 |
|
|
339 i = (*maxexp) + (*minexp); |
|
|
340 if (((*ibeta) == 2) && (i == 0)) (*maxexp) = (*maxexp) - 1; |
|
|
341 if (i > 20) (*maxexp) = (*maxexp) - 1; |
|
|
342 if (a != y) (*maxexp) = (*maxexp) - 2; |
|
|
343 (*xmax) = one - (*epsneg); |
|
|
344 tmp = (*xmax)*one; |
|
|
345 if (tmp != (*xmax)) (*xmax) = one - beta * (*epsneg); |
|
|
346 (*xmax) = (*xmax) / (beta * beta * beta * (*xmin)); |
|
|
347 i = (*maxexp) + (*minexp) + 3; |
|
|
348 if (i > 0) { |
|
|
349 |
|
|
350 for (j = 1; j<=i; j++ ) { |
|
|
351 if ((*ibeta) == 2) (*xmax) = (*xmax) + (*xmax); |
|
|
352 if ((*ibeta) != 2) (*xmax) = (*xmax) * beta; |
|
|
353 } |
|
|
354 |
|
|
355 } |
|
|
356 |
|
|
357 return; |
|
|
358 |
|
|
359 } |
|
|
360 |
|
|
361 typedef union |
|
|
362 { |
|
|
363 double d; |
|
|
364 int i[2]; |
|
|
365 } equiv; |
|
|
366 |
|
|
367 #ifdef DP |
|
|
368 int |
|
|
369 equiv_compare (equiv *std, equiv *v, int len) |
|
|
370 { |
|
|
371 int i; |
|
|
372 for (i = 0; i < len; i++) |
|
|
373 if (v[i].i[0] != std[i].i[0] || v[i].i[1] != std[i].i[1]) |
|
|
374 return 0; |
|
|
375 return 1; |
|
|
376 } |
|
|
377 #endif |
|
|
378 |
|
|
379 int |
|
|
380 main (void) |
|
|
381 { |
|
|
382 /* Works for 32 bit machines with 32 bit ints and 64 bit doubles */ |
|
|
383 |
|
|
384 int ibeta, iexp, irnd, it, machep, maxexp, minexp, negep, ngrd; |
|
|
385 REAL eps, epsneg, xmax, xmin; |
|
|
386 int i; |
|
|
387 equiv flt_params[4]; |
|
|
388 |
|
|
389 rmachar (&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp, |
|
|
390 &maxexp, &eps, &epsneg, &xmin, &xmax); |
|
|
391 |
|
|
392 flt_params[0].d = xmin; |
|
|
393 flt_params[1].d = xmax; |
|
|
394 flt_params[2].d = epsneg; |
|
|
395 flt_params[3].d = eps; |
|
|
396 |
|
|
397 #ifdef DP |
|
|
398 #define IS_MACH(v,nm,sm_1,sm_2,lrg_1,lrg_2,rt_1,rt_2,dv_1,dv_2) \ |
|
|
399 do \ |
|
|
400 { \ |
|
|
401 equiv v[4]; \ |
|
|
402 v[0].i[0] = (sm_1); v[0].i[1] = (sm_2); \ |
|
|
403 v[1].i[0] = (lrg_1); v[1].i[1] = (lrg_2); \ |
|
|
404 v[2].i[0] = (rt_1); v[2].i[1] = (rt_2); \ |
|
|
405 v[3].i[0] = (dv_1); v[3].i[1] = (dv_2); \ |
|
|
406 \ |
|
|
407 if (equiv_compare (v, flt_params, 4)) \ |
|
|
408 { \ |
|
|
409 printf ("%s\n", nm); \ |
|
|
410 return 0; \ |
|
|
411 } \ |
|
|
412 } \ |
|
|
413 while (0) |
|
|
414 |
|
|
415 IS_MACH (ieee_big_endian, "IEEE_BIG_ENDIAN", |
|
|
416 1048576, 0, |
|
|
417 2146435071, -1, |
|
|
418 1017118720, 0, |
|
|
419 1018167296, 0); |
|
|
420 /* 1070810131, 1352628735); */ |
|
|
421 |
|
|
422 IS_MACH (ieee_little_endian, "IEEE_LITTLE_ENDIAN", |
|
|
423 0, 1048576, |
|
|
424 -1, 2146435071, |
|
|
425 0, 1017118720, |
|
|
426 0, 1018167296); |
|
|
427 /* 1352628735, 1070810131); */ |
|
|
428 |
|
|
429 IS_MACH (vax_d_float, "VAX_D_FLOAT", |
|
|
430 128, 0, |
|
|
431 -32769, -1, |
|
|
432 9344, 0, |
|
|
433 9344, 0); |
|
|
434 /* 546979738, -805796613); */ |
|
|
435 |
|
|
436 IS_MACH (vax_g_float, "VAX_G_FLOAT", |
|
|
437 16, 0, |
|
|
438 -32769, -1, |
|
|
439 15552, 0, |
|
|
440 15552, 0); |
|
|
441 /* 1142112243, 2046775455); */ |
|
|
442 #else |
|
|
443 LOSE! LOSE! |
|
|
444 #endif |
|
|
445 |
|
|
446 printf ("UNRECOGNIZED_FLOATING_POINT_FORMAT\n"); |
|
|
447 return 1; |
|
|
448 } |
