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