140
141
142
143
144
145
146
147 LOGICAL IEEE1, RND
148 INTEGER BETA, T
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 LOGICAL FIRST, LIEEE1, LRND
193 INTEGER LBETA, LT
194 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
195
196
197 DOUBLE PRECISION DLAMC3
199
200
201 SAVE first, lieee1, lbeta, lrnd, lt
202
203
204 DATA first / .true. /
205
206
207
208 IF( first ) THEN
209 first = .false.
210 one = 1
211
212
213
214
215
216
217
218
219
220
221
222
223
224 a = 1
225 c = 1
226
227
228 10 CONTINUE
229 IF( c.EQ.one ) THEN
230 a = 2*a
233 GO TO 10
234 END IF
235
236
237
238
239
240
241
242 b = 1
244
245
246 20 CONTINUE
247 IF( c.EQ.a ) THEN
248 b = 2*b
250 GO TO 20
251 END IF
252
253
254
255
256
257
258
259 qtr = one / 4
260 savec = c
262 lbeta = c + qtr
263
264
265
266
267 b = lbeta
268 f =
dlamc3( b / 2, -b / 100 )
270 IF( c.EQ.a ) THEN
271 lrnd = .true.
272 ELSE
273 lrnd = .false.
274 END IF
275 f =
dlamc3( b / 2, b / 100 )
277 IF( ( lrnd ) .AND. ( c.EQ.a ) )
278 $ lrnd = .false.
279
280
281
282
283
284
285
287 t2 =
dlamc3( b / 2, savec )
288 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
289
290
291
292
293
294
295
296
297 lt = 0
298 a = 1
299 c = 1
300
301
302 30 CONTINUE
303 IF( c.EQ.one ) THEN
304 lt = lt + 1
305 a = a*lbeta
308 GO TO 30
309 END IF
310
311
312 END IF
313
314 beta = lbeta
315 t = lt
316 rnd = lrnd
317 ieee1 = lieee1
318 RETURN
319
320
321