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
|
double precision function nearfloat(x, dir)
*
* PURPOSE
* Compute the near (double) float from x in
* the direction dir
* dir >= 0 => toward +oo
* dir < 0 => toward -oo
*
* AUTHOR
* Bruno Pincon <Bruno.Pincon@iecn.u-nancy.fr>
*
* REMARK
* This code may be shorter if we assume that the
* radix base of floating point numbers is 2 : one
* could use the frexp C function to extract the
* mantissa and exponent part in place of dealing
* with a call to the log function with corrections
* to avoid possible floating point error...
*
implicit none
* PARAMETERS
double precision x, dir
* EXTERNAL FUNCTIONS
double precision dlamch
external dlamch
integer isanan
external isanan
* LOCAL VARIABLES
double precision sign_x, ep, xa, d, m
integer e, i, p
* STATIC VARIABLES
logical first, DENORM
double precision RMAX, RMIN, ULP, BASE, LNB, TINY
save RMAX, RMIN, ULP, BASE, LNB, TINY
save first, DENORM
data first /.true./
* TEXT
* got f.p. parameters used by the algorithm
if (first) then
RMAX = dlamch('o')
RMIN = dlamch('u')
BASE = dlamch('b')
p = int(dlamch('n'))
LNB = log(BASE)
* computation of 1 ulp : 1 ulp = base^(1-p)
* p = number of digits for the mantissa = dlamch('n')
ULP = BASE**(1 - p)
* query if denormalised numbers are used : if yes
* compute TINY the smallest denormalised number > 0 :
* TINY is also the increment between 2 neighbooring
* denormalised numbers
if (RMIN/BASE .ne. 0.d0) then
DENORM = .true.
TINY = RMIN
do i = 1, p-1
TINY = TINY / BASE
enddo
else
DENORM = .false.
endif
first = .false.
endif
d = sign(1.d0, dir)
sign_x = sign(1.d0, x)
xa = abs(x)
if (isanan(x) .eq. 1) then
* nan
nearfloat = x
elseif (xa .gt. RMAX) then
* +-inf
if (d*sign_x .lt. 0.d0) then
nearfloat = sign_x * RMAX
else
nearfloat = x
endif
elseif (xa .ge. RMIN) then
* usual case : x is a normalised floating point number
* 1/ got the exponent e and the exponent part ep = base^e
e = int( log(xa)/LNB )
ep = BASE**e
* in case of xa very near RMAX an error in e (of + 1)
* result in an overflow in ep
if ( ep .gt. RMAX ) then
e = e - 1
ep = BASE**e
endif
* also in case of xa very near RMIN and when denormalised
* number are not used, an error in e (of -1) results in a
* flush to 0 for ep.
if ( ep .eq. 0.d0 ) then
e = e + 1
ep = BASE**e
endif
* 2/ got the mantissa
m = xa/ep
* 3/ verify that 1 <= m < BASE
if ( m .lt. 1.d0 ) then
* multiply m by BASE
do while ( m .lt. 1.d0 )
m = m * BASE
ep = ep / BASE
enddo
elseif ( m .ge. BASE ) then
* divide m by BASE
do while ( m .ge. 1.d0 )
m = m / BASE
ep = ep * BASE
enddo
endif
* 4/ now compute the near float
if (d*sign_x .lt. 0.d0) then
* retrieve one ULP to m but there is a special case
if ( m .eq. 1.d0 ) then
* this is the special case : we must retrieve ULP / BASE
nearfloat = sign_x*( m - (ULP/BASE) )*ep
else
nearfloat = sign_x*( m - ULP )*ep
endif
else
nearfloat = sign_x*(m + ULP)*ep
endif
elseif ( x .eq. 0.d0 ) then
* case x = 0 nearfloat depends if denormalised numbers are used
if (DENORM) then
nearfloat = d*TINY
else
nearfloat = d*RMIN
endif
else
* x is a denormalised number
nearfloat = x + d*TINY
endif
end
|
