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
|
*
************************************************************************
*
SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
*
* -- LAPACK auxiliary routine (version 2.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* October 31, 1992
*
* .. Scalar Arguments ..
LOGICAL IEEE1, RND
INTEGER BETA, T
* ..
*
* Purpose
* =======
*
* DLAMC1 determines the machine parameters given by BETA, T, RND, and
* IEEE1.
*
* Arguments
* =========
*
* BETA (output) INTEGER
* The base of the machine.
*
* T (output) INTEGER
* The number of ( BETA ) digits in the mantissa.
*
* RND (output) LOGICAL
* Specifies whether proper rounding ( RND = .TRUE. ) or
* chopping ( RND = .FALSE. ) occurs in addition. This may not
* be a reliable guide to the way in which the machine performs
* its arithmetic.
*
* IEEE1 (output) LOGICAL
* Specifies whether rounding appears to be done in the IEEE
* 'round to nearest' style.
*
* Further Details
* ===============
*
* The routine is based on the routine ENVRON by Malcolm and
* incorporates suggestions by Gentleman and Marovich. See
*
* Malcolm M. A. (1972) Algorithms to reveal properties of
* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
*
* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
* that reveal properties of floating point arithmetic units.
* Comms. of the ACM, 17, 276-277.
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL FIRST, LIEEE1, LRND
INTEGER LBETA, LT
DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMC3
EXTERNAL DLAMC3
* ..
* .. Save statement ..
SAVE FIRST, LIEEE1, LBETA, LRND, LT
* ..
* .. Data statements ..
DATA FIRST / .TRUE. /
* ..
* .. Executable Statements ..
*
IF( FIRST ) THEN
FIRST = .FALSE.
ONE = 1
*
* LBETA, LIEEE1, LT and LRND are the local values of BETA,
* IEEE1, T and RND.
*
* Throughout this routine we use the function DLAMC3 to ensure
* that relevant values are stored and not held in registers, or
* are not affected by optimizers.
*
* Compute a = 2.0**m with the smallest positive integer m such
* that
*
* fl( a + 1.0 ) = a.
*
A = 1
C = 1
*
*+ WHILE( C.EQ.ONE )LOOP
10 CONTINUE
IF( C.EQ.ONE ) THEN
A = 2*A
C = DLAMC3( A, ONE )
C = DLAMC3( C, -A )
GO TO 10
END IF
*+ END WHILE
*
* Now compute b = 2.0**m with the smallest positive integer m
* such that
*
* fl( a + b ) .gt. a.
*
B = 1
C = DLAMC3( A, B )
*
*+ WHILE( C.EQ.A )LOOP
20 CONTINUE
IF( C.EQ.A ) THEN
B = 2*B
C = DLAMC3( A, B )
GO TO 20
END IF
*+ END WHILE
*
* Now compute the base. a and c are neighbouring floating point
* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
* their difference is beta. Adding 0.25 to c is to ensure that it
* is truncated to beta and not ( beta - 1 ).
*
QTR = ONE / 4
SAVEC = C
C = DLAMC3( C, -A )
LBETA = C + QTR
*
* Now determine whether rounding or chopping occurs, by adding a
* bit less than beta/2 and a bit more than beta/2 to a.
*
B = LBETA
F = DLAMC3( B / 2, -B / 100 )
C = DLAMC3( F, A )
IF( C.EQ.A ) THEN
LRND = .TRUE.
ELSE
LRND = .FALSE.
END IF
F = DLAMC3( B / 2, B / 100 )
C = DLAMC3( F, A )
IF( ( LRND ) .AND. ( C.EQ.A ) )
$ LRND = .FALSE.
*
* Try and decide whether rounding is done in the IEEE 'round to
* nearest' style. B/2 is half a unit in the last place of the two
* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
* zero, and SAVEC is odd. Thus adding B/2 to A should not change
* A, but adding B/2 to SAVEC should change SAVEC.
*
T1 = DLAMC3( B / 2, A )
T2 = DLAMC3( B / 2, SAVEC )
LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
*
* Now find the mantissa, t. It should be the integer part of
* log to the base beta of a, however it is safer to determine t
* by powering. So we find t as the smallest positive integer for
* which
*
* fl( beta**t + 1.0 ) = 1.0.
*
LT = 0
A = 1
C = 1
*
*+ WHILE( C.EQ.ONE )LOOP
30 CONTINUE
IF( C.EQ.ONE ) THEN
LT = LT + 1
A = A*LBETA
C = DLAMC3( A, ONE )
C = DLAMC3( C, -A )
GO TO 30
END IF
*+ END WHILE
*
END IF
*
BETA = LBETA
T = LT
RND = LRND
IEEE1 = LIEEE1
RETURN
*
* End of DLAMC1
*
END
|
