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