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