|
173
|
1 !===-- module/ieee_arithmetic.f90 ------------------------------------------===!
|
|
|
2 !
|
|
|
3 ! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
|
|
|
4 ! See https://llvm.org/LICENSE.txt for license information.
|
|
|
5 ! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
|
|
|
6 !
|
|
|
7 !===------------------------------------------------------------------------===!
|
|
|
8
|
|
|
9 ! See Fortran 2018, clause 17.2
|
|
|
10 module ieee_arithmetic
|
|
|
11
|
|
221
|
12 use __Fortran_builtins, only: &
|
|
|
13 ieee_is_nan => __builtin_ieee_is_nan, &
|
|
|
14 ieee_selected_real_kind => __builtin_ieee_selected_real_kind, &
|
|
|
15 ieee_support_datatype => __builtin_ieee_support_datatype, &
|
|
|
16 ieee_support_denormal => __builtin_ieee_support_denormal, &
|
|
|
17 ieee_support_divide => __builtin_ieee_support_divide, &
|
|
|
18 ieee_support_inf => __builtin_ieee_support_inf, &
|
|
|
19 ieee_support_io => __builtin_ieee_support_io, &
|
|
|
20 ieee_support_nan => __builtin_ieee_support_nan, &
|
|
|
21 ieee_support_sqrt => __builtin_ieee_support_sqrt, &
|
|
|
22 ieee_support_standard => __builtin_ieee_support_standard, &
|
|
|
23 ieee_support_subnormal => __builtin_ieee_support_subnormal, &
|
|
|
24 ieee_support_underflow_control => __builtin_ieee_support_underflow_control
|
|
|
25
|
|
|
26 implicit none
|
|
|
27
|
|
173
|
28 type :: ieee_class_type
|
|
|
29 private
|
|
|
30 integer(kind=1) :: which = 0
|
|
|
31 end type ieee_class_type
|
|
|
32
|
|
|
33 type(ieee_class_type), parameter :: &
|
|
|
34 ieee_signaling_nan = ieee_class_type(1), &
|
|
|
35 ieee_quiet_nan = ieee_class_type(2), &
|
|
|
36 ieee_negative_inf = ieee_class_type(3), &
|
|
|
37 ieee_negative_normal = ieee_class_type(4), &
|
|
|
38 ieee_negative_denormal = ieee_class_type(5), &
|
|
|
39 ieee_negative_zero = ieee_class_type(6), &
|
|
|
40 ieee_positive_zero = ieee_class_type(7), &
|
|
|
41 ieee_positive_subnormal = ieee_class_type(8), &
|
|
|
42 ieee_positive_normal = ieee_class_type(9), &
|
|
|
43 ieee_positive_inf = ieee_class_type(10), &
|
|
|
44 ieee_other_value = ieee_class_type(11)
|
|
|
45
|
|
|
46 type(ieee_class_type), parameter :: &
|
|
|
47 ieee_negative_subnormal = ieee_negative_denormal, &
|
|
|
48 ieee_positive_denormal = ieee_negative_subnormal
|
|
|
49
|
|
|
50 type :: ieee_round_type
|
|
|
51 private
|
|
|
52 integer(kind=1) :: mode = 0
|
|
|
53 end type ieee_round_type
|
|
|
54
|
|
|
55 type(ieee_round_type), parameter :: &
|
|
|
56 ieee_nearest = ieee_round_type(1), &
|
|
|
57 ieee_to_zero = ieee_round_type(2), &
|
|
|
58 ieee_up = ieee_round_type(3), &
|
|
|
59 ieee_down = ieee_round_type(4), &
|
|
|
60 ieee_away = ieee_round_type(5), &
|
|
|
61 ieee_other = ieee_round_type(6)
|
|
|
62
|
|
|
63 interface operator(==)
|
|
|
64 module procedure class_eq
|
|
|
65 module procedure round_eq
|
|
|
66 end interface operator(==)
|
|
|
67 interface operator(/=)
|
|
|
68 module procedure class_ne
|
|
|
69 module procedure round_ne
|
|
|
70 end interface operator(/=)
|
|
221
|
71 private :: class_eq, class_ne, round_eq, round_ne
|
|
173
|
72
|
|
|
73 ! See Fortran 2018, 17.10 & 17.11
|
|
221
|
74 generic :: ieee_class => ieee_class_a2, ieee_class_a3, ieee_class_a4, ieee_class_a8, ieee_class_a10, ieee_class_a16
|
|
|
75 private :: ieee_class_a2, ieee_class_a3, ieee_class_a4, ieee_class_a8, ieee_class_a10, ieee_class_a16
|
|
|
76
|
|
|
77 generic :: ieee_copy_sign => ieee_copy_sign_a2, ieee_copy_sign_a3, ieee_copy_sign_a4, ieee_copy_sign_a8, ieee_copy_sign_a10, ieee_copy_sign_a16
|
|
|
78 private :: ieee_copy_sign_a2, ieee_copy_sign_a3, ieee_copy_sign_a4, ieee_copy_sign_a8, ieee_copy_sign_a10, ieee_copy_sign_a16
|
|
|
79
|
|
|
80 generic :: ieee_is_finite => ieee_is_finite_a2, ieee_is_finite_a3, ieee_is_finite_a4, ieee_is_finite_a8, ieee_is_finite_a10, ieee_is_finite_a16
|
|
|
81 private :: ieee_is_finite_a2, ieee_is_finite_a3, ieee_is_finite_a4, ieee_is_finite_a8, ieee_is_finite_a10, ieee_is_finite_a16
|
|
173
|
82
|
|
221
|
83 generic :: ieee_rem => &
|
|
|
84 ieee_rem_a2_a2, ieee_rem_a2_a3, ieee_rem_a2_a4, ieee_rem_a2_a8, ieee_rem_a2_a10, ieee_rem_a2_a16, &
|
|
|
85 ieee_rem_a3_a2, ieee_rem_a3_a3, ieee_rem_a3_a4, ieee_rem_a3_a8, ieee_rem_a3_a10, ieee_rem_a3_a16, &
|
|
|
86 ieee_rem_a4_a2, ieee_rem_a4_a3, ieee_rem_a4_a4, ieee_rem_a4_a8, ieee_rem_a4_a10, ieee_rem_a4_a16, &
|
|
|
87 ieee_rem_a8_a2, ieee_rem_a8_a3, ieee_rem_a8_a4, ieee_rem_a8_a8, ieee_rem_a8_a10, ieee_rem_a8_a16, &
|
|
|
88 ieee_rem_a10_a2, ieee_rem_a10_a3, ieee_rem_a10_a4, ieee_rem_a10_a8, ieee_rem_a10_a10, ieee_rem_a10_a16, &
|
|
|
89 ieee_rem_a16_a2, ieee_rem_a16_a3, ieee_rem_a16_a4, ieee_rem_a16_a8, ieee_rem_a16_a10, ieee_rem_a16_a16
|
|
|
90 private :: &
|
|
|
91 ieee_rem_a2_a2, ieee_rem_a2_a3, ieee_rem_a2_a4, ieee_rem_a2_a8, ieee_rem_a2_a10, ieee_rem_a2_a16, &
|
|
|
92 ieee_rem_a3_a2, ieee_rem_a3_a3, ieee_rem_a3_a4, ieee_rem_a3_a8, ieee_rem_a3_a10, ieee_rem_a3_a16, &
|
|
|
93 ieee_rem_a4_a2, ieee_rem_a4_a3, ieee_rem_a4_a4, ieee_rem_a4_a8, ieee_rem_a4_a10, ieee_rem_a4_a16, &
|
|
|
94 ieee_rem_a8_a2, ieee_rem_a8_a3, ieee_rem_a8_a4, ieee_rem_a8_a8, ieee_rem_a8_a10, ieee_rem_a8_a16, &
|
|
|
95 ieee_rem_a10_a2, ieee_rem_a10_a3, ieee_rem_a10_a4, ieee_rem_a10_a8, ieee_rem_a10_a10, ieee_rem_a10_a16, &
|
|
|
96 ieee_rem_a16_a2, ieee_rem_a16_a3, ieee_rem_a16_a4, ieee_rem_a16_a8, ieee_rem_a16_a10, ieee_rem_a16_a16
|
|
|
97
|
|
|
98 generic :: ieee_support_rounding => ieee_support_rounding_, &
|
|
|
99 ieee_support_rounding_2, ieee_support_rounding_3, &
|
|
|
100 ieee_support_rounding_4, ieee_support_rounding_8, &
|
|
|
101 ieee_support_rounding_10, ieee_support_rounding_16
|
|
|
102 private :: ieee_support_rounding_, &
|
|
|
103 ieee_support_rounding_2, ieee_support_rounding_3, &
|
|
|
104 ieee_support_rounding_4, ieee_support_rounding_8, &
|
|
|
105 ieee_support_rounding_10, ieee_support_rounding_16
|
|
173
|
106
|
|
|
107 ! TODO: more interfaces (_fma, &c.)
|
|
|
108
|
|
|
109 private :: classify
|
|
|
110
|
|
|
111 contains
|
|
|
112
|
|
|
113 elemental logical function class_eq(x,y)
|
|
|
114 type(ieee_class_type), intent(in) :: x, y
|
|
|
115 class_eq = x%which == y%which
|
|
|
116 end function class_eq
|
|
|
117
|
|
|
118 elemental logical function class_ne(x,y)
|
|
|
119 type(ieee_class_type), intent(in) :: x, y
|
|
|
120 class_ne = x%which /= y%which
|
|
|
121 end function class_ne
|
|
|
122
|
|
|
123 elemental logical function round_eq(x,y)
|
|
|
124 type(ieee_round_type), intent(in) :: x, y
|
|
|
125 round_eq = x%mode == y%mode
|
|
|
126 end function round_eq
|
|
|
127
|
|
|
128 elemental logical function round_ne(x,y)
|
|
|
129 type(ieee_round_type), intent(in) :: x, y
|
|
|
130 round_ne = x%mode /= y%mode
|
|
|
131 end function round_ne
|
|
|
132
|
|
|
133 elemental type(ieee_class_type) function classify( &
|
|
|
134 expo,maxExpo,negative,significandNZ,quietBit)
|
|
|
135 integer, intent(in) :: expo, maxExpo
|
|
|
136 logical, intent(in) :: negative, significandNZ, quietBit
|
|
|
137 if (expo == 0) then
|
|
|
138 if (significandNZ) then
|
|
|
139 if (negative) then
|
|
|
140 classify = ieee_negative_denormal
|
|
|
141 else
|
|
|
142 classify = ieee_positive_denormal
|
|
|
143 end if
|
|
|
144 else
|
|
|
145 if (negative) then
|
|
|
146 classify = ieee_negative_zero
|
|
|
147 else
|
|
|
148 classify = ieee_positive_zero
|
|
|
149 end if
|
|
|
150 end if
|
|
|
151 else if (expo == maxExpo) then
|
|
|
152 if (significandNZ) then
|
|
|
153 if (quietBit) then
|
|
|
154 classify = ieee_quiet_nan
|
|
|
155 else
|
|
|
156 classify = ieee_signaling_nan
|
|
|
157 end if
|
|
|
158 else
|
|
|
159 if (negative) then
|
|
|
160 classify = ieee_negative_inf
|
|
|
161 else
|
|
|
162 classify = ieee_positive_inf
|
|
|
163 end if
|
|
|
164 end if
|
|
|
165 else
|
|
|
166 if (negative) then
|
|
|
167 classify = ieee_negative_normal
|
|
|
168 else
|
|
|
169 classify = ieee_positive_normal
|
|
|
170 end if
|
|
|
171 end if
|
|
|
172 end function classify
|
|
|
173
|
|
|
174 #define _CLASSIFY(RKIND,IKIND,TOTALBITS,PREC,IMPLICIT) \
|
|
|
175 type(ieee_class_type) elemental function ieee_class_a##RKIND(x); \
|
|
|
176 real(kind=RKIND), intent(in) :: x; \
|
|
|
177 integer(kind=IKIND) :: raw; \
|
|
|
178 integer, parameter :: significand = PREC - IMPLICIT; \
|
|
|
179 integer, parameter :: exponentBits = TOTALBITS - 1 - significand; \
|
|
|
180 integer, parameter :: maxExpo = shiftl(1, exponentBits) - 1; \
|
|
|
181 integer :: exponent, sign; \
|
|
|
182 logical :: negative, nzSignificand, quiet; \
|
|
|
183 raw = transfer(x, raw); \
|
|
|
184 exponent = ibits(raw, significand, exponentBits); \
|
|
|
185 negative = btest(raw, TOTALBITS - 1); \
|
|
|
186 nzSignificand = ibits(raw, 0, significand) /= 0; \
|
|
|
187 quiet = btest(raw, significand - 1); \
|
|
|
188 ieee_class_a##RKIND = classify(exponent, maxExpo, negative, nzSignificand, quiet); \
|
|
|
189 end function ieee_class_a##RKIND
|
|
|
190 _CLASSIFY(2,2,16,11,1)
|
|
|
191 _CLASSIFY(3,2,16,8,1)
|
|
|
192 _CLASSIFY(4,4,32,24,1)
|
|
|
193 _CLASSIFY(8,8,64,53,1)
|
|
|
194 _CLASSIFY(10,16,80,64,0)
|
|
|
195 _CLASSIFY(16,16,128,112,1)
|
|
|
196 #undef _CLASSIFY
|
|
|
197
|
|
|
198 ! TODO: This might need to be an actual Operation instead
|
|
|
199 #define _COPYSIGN(RKIND,IKIND,BITS) \
|
|
|
200 real(kind=RKIND) elemental function ieee_copy_sign_a##RKIND(x,y); \
|
|
|
201 real(kind=RKIND), intent(in) :: x, y; \
|
|
|
202 integer(kind=IKIND) :: xbits, ybits; \
|
|
|
203 xbits = transfer(x, xbits); \
|
|
|
204 ybits = transfer(y, ybits); \
|
|
|
205 xbits = ior(ibclr(xbits, BITS-1), iand(ybits, shiftl(1_##IKIND, BITS-1))); \
|
|
|
206 ieee_copy_sign_a##RKIND = transfer(xbits, x); \
|
|
|
207 end function ieee_copy_sign_a##RKIND
|
|
|
208 _COPYSIGN(2,2,16)
|
|
|
209 _COPYSIGN(3,2,16)
|
|
|
210 _COPYSIGN(4,4,32)
|
|
|
211 _COPYSIGN(8,8,64)
|
|
|
212 _COPYSIGN(10,16,80)
|
|
|
213 _COPYSIGN(16,16,128)
|
|
|
214 #undef _COPYSIGN
|
|
|
215
|
|
221
|
216 #define _IS_FINITE(KIND) \
|
|
|
217 elemental function ieee_is_finite_a##KIND(x) result(res); \
|
|
|
218 real(kind=KIND), intent(in) :: x; \
|
|
|
219 logical :: res; \
|
|
|
220 type(ieee_class_type) :: classification; \
|
|
|
221 classification = ieee_class(x); \
|
|
|
222 res = classification == ieee_negative_zero .or. classification == ieee_positive_zero \
|
|
|
223 .or. classification == ieee_negative_denormal .or. classification == ieee_positive_denormal \
|
|
|
224 .or. classification == ieee_negative_normal .or. classification == ieee_positive_normal; \
|
|
|
225 end function
|
|
|
226 _IS_FINITE(2)
|
|
|
227 _IS_FINITE(3)
|
|
|
228 _IS_FINITE(4)
|
|
|
229 _IS_FINITE(8)
|
|
|
230 _IS_FINITE(10)
|
|
|
231 _IS_FINITE(16)
|
|
|
232 #undef _IS_FINITE
|
|
|
233
|
|
|
234 ! TODO: handle edge cases from 17.11.31
|
|
|
235 #define _REM(XKIND,YKIND) \
|
|
|
236 elemental function ieee_rem_a##XKIND##_a##YKIND(x, y) result(res); \
|
|
|
237 real(kind=XKIND), intent(in) :: x; \
|
|
|
238 real(kind=YKIND), intent(in) :: y; \
|
|
|
239 integer, parameter :: rkind = max(XKIND, YKIND); \
|
|
|
240 real(kind=rkind) :: res, tmp; \
|
|
|
241 tmp = anint(real(x, kind=rkind) / y); \
|
|
|
242 res = x - y * tmp; \
|
|
|
243 end function
|
|
|
244 _REM(2,2)
|
|
|
245 _REM(2,3)
|
|
|
246 _REM(2,4)
|
|
|
247 _REM(2,8)
|
|
|
248 _REM(2,10)
|
|
|
249 _REM(2,16)
|
|
|
250 _REM(3,2)
|
|
|
251 _REM(3,3)
|
|
|
252 _REM(3,4)
|
|
|
253 _REM(3,8)
|
|
|
254 _REM(3,10)
|
|
|
255 _REM(3,16)
|
|
|
256 _REM(4,2)
|
|
|
257 _REM(4,3)
|
|
|
258 _REM(4,4)
|
|
|
259 _REM(4,8)
|
|
|
260 _REM(4,10)
|
|
|
261 _REM(4,16)
|
|
|
262 _REM(8,2)
|
|
|
263 _REM(8,3)
|
|
|
264 _REM(8,4)
|
|
|
265 _REM(8,8)
|
|
|
266 _REM(8,10)
|
|
|
267 _REM(8,16)
|
|
|
268 _REM(10,2)
|
|
|
269 _REM(10,3)
|
|
|
270 _REM(10,4)
|
|
|
271 _REM(10,8)
|
|
|
272 _REM(10,10)
|
|
|
273 _REM(10,16)
|
|
|
274 _REM(16,2)
|
|
|
275 _REM(16,3)
|
|
|
276 _REM(16,4)
|
|
|
277 _REM(16,8)
|
|
|
278 _REM(16,10)
|
|
|
279 _REM(16,16)
|
|
|
280 #undef _REM
|
|
|
281
|
|
|
282 pure logical function ieee_support_rounding_(round_type)
|
|
|
283 type(ieee_round_type), intent(in) :: round_type
|
|
|
284 ieee_support_rounding_ = .true.
|
|
|
285 end function
|
|
|
286 pure logical function ieee_support_rounding_2(round_type,x)
|
|
|
287 type(ieee_round_type), intent(in) :: round_type
|
|
|
288 real(kind=2), intent(in) :: x
|
|
|
289 ieee_support_rounding_2 = .true.
|
|
|
290 end function
|
|
|
291 pure logical function ieee_support_rounding_3(round_type,x)
|
|
|
292 type(ieee_round_type), intent(in) :: round_type
|
|
|
293 real(kind=3), intent(in) :: x
|
|
|
294 ieee_support_rounding_3 = .true.
|
|
|
295 end function
|
|
|
296 pure logical function ieee_support_rounding_4(round_type,x)
|
|
|
297 type(ieee_round_type), intent(in) :: round_type
|
|
|
298 real(kind=4), intent(in) :: x
|
|
|
299 ieee_support_rounding_4 = .true.
|
|
|
300 end function
|
|
|
301 pure logical function ieee_support_rounding_8(round_type,x)
|
|
|
302 type(ieee_round_type), intent(in) :: round_type
|
|
|
303 real(kind=8), intent(in) :: x
|
|
|
304 ieee_support_rounding_8 = .true.
|
|
|
305 end function
|
|
|
306 pure logical function ieee_support_rounding_10(round_type,x)
|
|
|
307 type(ieee_round_type), intent(in) :: round_type
|
|
|
308 real(kind=10), intent(in) :: x
|
|
|
309 ieee_support_rounding_10 = .true.
|
|
|
310 end function
|
|
|
311 pure logical function ieee_support_rounding_16(round_type,x)
|
|
|
312 type(ieee_round_type), intent(in) :: round_type
|
|
|
313 real(kind=16), intent(in) :: x
|
|
|
314 ieee_support_rounding_16 = .true.
|
|
|
315 end function
|
|
|
316
|
|
173
|
317 end module ieee_arithmetic
|
