1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
|
cc
c
cc
Double precision function dbl_eps ()
integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
double precision eps,epsneg,xmin,xmax
call machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,
$ eps,epsneg,xmin,xmax)
dbl_eps = eps / 2.0d0
return
end
cc
c
cc
Double precision function dbl_epsneg ()
integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
double precision eps,epsneg,xmin,xmax
call machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,
$ eps,epsneg,xmin,xmax)
dbl_epsneg = epsneg
return
end
cc
c
cc
Double precision function dbl_max ()
integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
double precision deps,epsneg,xmin,xmax
call machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,
$ deps,epsneg,xmin,xmax)
dbl_max = xmax
return
end
cc
c
cc
Integer function dbl_max_exp ()
integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
double precision deps,epsneg,xmin,xmax
call machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,
$ deps,epsneg,xmin,xmax)
dbl_max_exp = Int(Dlog10(xmax))
return
end
cc
c
cc
Double precision function dbl_min ()
integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
double precision deps,epsneg,xmin,xmax
call machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,
$ deps,epsneg,xmin,xmax)
dbl_min = xmin
return
end
cc
c
cc
Integer function dbl_min_exp ()
integer ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp
double precision deps,epsneg,xmin,xmax
call machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,maxexp,
$ deps,epsneg,xmin,xmax)
dbl_min_exp = Int(Dlog10(xmin))
return
end
C ALGORITHM 665, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 14, NO. 4, PP. 303-311.
SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
1 MAXEXP,EPS,EPSNEG,XMIN,XMAX)
C-----------------------------------------------------------------------
C This Fortran 77 subroutine is intended to determine the parameters
C of the floating-point arithmetic system specified below. The
C determination of the first three uses an extension of an algorithm
C due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
C but not all, of the improvements suggested by M. Gentleman and S.
C Marovich, CACM 17 (1974), pp. 276-277. An earlier version of this
C program was published in the book Software Manual for the
C Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
C Englewood Cliffs, NJ, 1980.
C
C The program as given here must be modified before compiling. If
C a single (double) precision version is desired, change all
C occurrences of CS (CD) in columns 1 and 2 to blanks.
C
C Parameter values reported are as follows:
C
C IBETA - the radix for the floating-point representation
C IT - the number of base IBETA digits in the floating-point
C significand
C IRND - 0 if floating-point addition chops
C 1 if floating-point addition rounds, but not in the
C IEEE style
C 2 if floating-point addition rounds in the IEEE style
C 3 if floating-point addition chops, and there is
C partial underflow
C 4 if floating-point addition rounds, but not in the
C IEEE style, and there is partial underflow
C 5 if floating-point addition rounds in the IEEE style,
C and there is partial underflow
C NGRD - the number of guard digits for multiplication with
C truncating arithmetic. It is
C 0 if floating-point arithmetic rounds, or if it
C truncates and only IT base IBETA digits
C participate in the post-normalization shift of the
C floating-point significand in multiplication;
C 1 if floating-point arithmetic truncates and more
C than IT base IBETA digits participate in the
C post-normalization shift of the floating-point
C significand in multiplication.
C MACHEP - the largest negative integer such that
C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
C MACHEP is bounded below by -(IT+3)
C NEGEPS - the largest negative integer such that
C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
C NEGEPS is bounded below by -(IT+3)
C IEXP - the number of bits (decimal places if IBETA = 10)
C reserved for the representation of the exponent
C (including the bias or sign) of a floating-point
C number
C MINEXP - the largest in magnitude negative integer such that
C FLOAT(IBETA)**MINEXP is positive and normalized
C MAXEXP - the smallest positive power of BETA that overflows
C EPS - the smallest positive floating-point number such
C that 1.0+EPS .NE. 1.0. In particular, if either
C IBETA = 2 or IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
C Otherwise, EPS = (FLOAT(IBETA)**MACHEP)/2
C EPSNEG - A small positive floating-point number such that
C 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
C or IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
C Otherwise, EPSNEG = (IBETA**NEGEPS)/2. Because
C NEGEPS is bounded below by -(IT+3), EPSNEG may not
C be the smallest number that can alter 1.0 by
C subtraction.
C XMIN - the smallest non-vanishing normalized floating-point
C power of the radix, i.e., XMIN = FLOAT(IBETA)**MINEXP
C XMAX - the largest finite floating-point number. In
C particular XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
C Note - on some machines XMAX will be only the
C second, or perhaps third, largest number, being
C too small by 1 or 2 units in the last digit of
C the significand.
C
C Latest revision - April 20, 1987
C
C Author - W. J. Cody
C Argonne National Laboratory
C
C-----------------------------------------------------------------------
INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP,
1 MINEXP,MX,NEGEP,NGRD,NXRES
CS REAL A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA,
CS 1 TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
DOUBLE PRECISION A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,
1 T,TEMP,TEMPA,TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
C-----------------------------------------------------------------------
CS CONV(I) = REAL(I)
CONV(I) = DBLE(I)
ONE = CONV(1)
TWO = ONE + ONE
ZERO = ONE - ONE
C-----------------------------------------------------------------------
C Determine IBETA, BETA ala Malcolm.
C-----------------------------------------------------------------------
A = ONE
10 A = A + A
TEMP = A+ONE
TEMP1 = TEMP-A
IF (TEMP1-ONE .EQ. ZERO) GO TO 10
B = ONE
20 B = B + B
TEMP = A+B
ITEMP = INT(TEMP-A)
IF (ITEMP .EQ. 0) GO TO 20
IBETA = ITEMP
BETA = CONV(IBETA)
C-----------------------------------------------------------------------
C Determine IT, IRND.
C-----------------------------------------------------------------------
IT = 0
B = ONE
100 IT = IT + 1
B = B * BETA
TEMP = B+ONE
TEMP1 = TEMP-B
IF (TEMP1-ONE .EQ. ZERO) GO TO 100
IRND = 0
BETAH = BETA / TWO
TEMP = A+BETAH
IF (TEMP-A .NE. ZERO) IRND = 1
TEMPA = A + BETA
TEMP = TEMPA+BETAH
IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2
C-----------------------------------------------------------------------
C Determine NEGEP, EPSNEG.
C-----------------------------------------------------------------------
NEGEP = IT + 3
BETAIN = ONE / BETA
A = ONE
DO 200 I = 1, NEGEP
A = A * BETAIN
200 CONTINUE
B = A
210 TEMP = ONE-A
IF (TEMP-ONE .NE. ZERO) GO TO 220
A = A * BETA
NEGEP = NEGEP - 1
GO TO 210
220 NEGEP = -NEGEP
EPSNEG = A
IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 300
A = (A*(ONE+A)) / TWO
TEMP = ONE-A
IF (TEMP-ONE .NE. ZERO) EPSNEG = A
C-----------------------------------------------------------------------
C Determine MACHEP, EPS.
C-----------------------------------------------------------------------
300 MACHEP = -IT - 3
A = B
310 TEMP = ONE+A
IF (TEMP-ONE .NE. ZERO) GO TO 320
A = A * BETA
MACHEP = MACHEP + 1
GO TO 310
320 EPS = A
TEMP = TEMPA+BETA*(ONE+EPS)
IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 350
A = (A*(ONE+A)) / TWO
TEMP = ONE+A
IF (TEMP-ONE .NE. ZERO) EPS = A
C-----------------------------------------------------------------------
C Determine NGRD.
C-----------------------------------------------------------------------
350 NGRD = 0
TEMP = ONE+EPS
IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1
C-----------------------------------------------------------------------
C Determine IEXP, MINEXP, XMIN.
C
C Loop to determine largest I and K = 2**I such that
C (1/BETA) ** (2**(I))
C does not underflow.
C Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------------
I = 0
K = 1
Z = BETAIN
T = ONE + EPS
NXRES = 0
400 Y = Z
Z = Y * Y
C-----------------------------------------------------------------------
C Check for underflow here.
C-----------------------------------------------------------------------
A = Z * ONE
TEMP = Z * T
IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410
TEMP1 = TEMP * BETAIN
IF (TEMP1*BETA .EQ. Z) GO TO 410
I = I + 1
K = K + K
GO TO 400
410 IF (IBETA .EQ. 10) GO TO 420
IEXP = I + 1
MX = K + K
GO TO 450
C-----------------------------------------------------------------------
C This segment is for decimal machines only.
C-----------------------------------------------------------------------
420 IEXP = 2
IZ = IBETA
430 IF (K .LT. IZ) GO TO 440
IZ = IZ * IBETA
IEXP = IEXP + 1
GO TO 430
440 MX = IZ + IZ - 1
C-----------------------------------------------------------------------
C Loop to determine MINEXP, XMIN.
C Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------------
450 XMIN = Y
C-----------------------------------------------------------------------
C Check for underflow here.
C-----------------------------------------------------------------------
c IF (((Y*ONE*BETAIN+Y*ONE*BETAIN) .EQ. ZERO)
c 1 .OR. (ABS(Y*BETAIN) .GE. XMIN)) GO TO 460
Y = Y * BETAIN
A = Y * ONE
IF (((A+A) .EQ. ZERO)
1 .OR. (ABS(A) .GE. XMIN)) GO TO 460
TEMP = Y * T
K = K + 1
TEMP1 = TEMP * BETAIN
IF (TEMP1*BETA .NE. Y) GO TO 450
NXRES = 3
XMIN = Y
460 MINEXP = -K
C-----------------------------------------------------------------------
C Determine MAXEXP, XMAX.
C-----------------------------------------------------------------------
IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500
MX = MX + MX
IEXP = IEXP + 1
500 MAXEXP = MX + MINEXP
C-----------------------------------------------------------------
C Adjust IRND to reflect partial underflow.
C-----------------------------------------------------------------
IRND = IRND + NXRES
C-----------------------------------------------------------------
C Adjust for IEEE-style machines.
C-----------------------------------------------------------------
IF ((IRND .EQ. 2) .OR. (IRND .EQ. 5)) MAXEXP = MAXEXP - 2
C-----------------------------------------------------------------
C Adjust for non-IEEE machines with partial underflow.
C-----------------------------------------------------------------
IF ((IRND .EQ. 3) .OR. (IRND .EQ. 4)) MAXEXP = MAXEXP - IT
C-----------------------------------------------------------------
C Adjust for machines with implicit leading bit in binary
C significand, and machines with radix point at extreme
C right of significand.
C-----------------------------------------------------------------
I = MAXEXP + MINEXP
IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1
IF (I .GT. 20) MAXEXP = MAXEXP - 1
IF (A .NE. Y) MAXEXP = MAXEXP - 2
XMAX = ONE - EPSNEG
IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
XMAX = XMAX / (BETA * BETA * BETA * XMIN)
I = MAXEXP + MINEXP + 3
IF (I .LE. 0) GO TO 520
DO 510 J = 1, I
IF (IBETA .EQ. 2) XMAX = XMAX + XMAX
IF (IBETA .NE. 2) XMAX = XMAX * BETA
510 CONTINUE
520 RETURN
C---------- LAST CARD OF MACHAR ----------
END
|
