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 type :: ieee_class_type 13 private 14 integer(kind=1) :: which = 0 15 end type ieee_class_type 16 17 type(ieee_class_type), parameter :: & 18 ieee_signaling_nan = ieee_class_type(1), & 19 ieee_quiet_nan = ieee_class_type(2), & 20 ieee_negative_inf = ieee_class_type(3), & 21 ieee_negative_normal = ieee_class_type(4), & 22 ieee_negative_denormal = ieee_class_type(5), & 23 ieee_negative_zero = ieee_class_type(6), & 24 ieee_positive_zero = ieee_class_type(7), & 25 ieee_positive_subnormal = ieee_class_type(8), & 26 ieee_positive_normal = ieee_class_type(9), & 27 ieee_positive_inf = ieee_class_type(10), & 28 ieee_other_value = ieee_class_type(11) 29 30 type(ieee_class_type), parameter :: & 31 ieee_negative_subnormal = ieee_negative_denormal, & 32 ieee_positive_denormal = ieee_negative_subnormal 33 34 type :: ieee_round_type 35 private 36 integer(kind=1) :: mode = 0 37 end type ieee_round_type 38 39 type(ieee_round_type), parameter :: & 40 ieee_nearest = ieee_round_type(1), & 41 ieee_to_zero = ieee_round_type(2), & 42 ieee_up = ieee_round_type(3), & 43 ieee_down = ieee_round_type(4), & 44 ieee_away = ieee_round_type(5), & 45 ieee_other = ieee_round_type(6) 46 47 interface operator(==) 48 module procedure class_eq 49 module procedure round_eq 50 end interface operator(==) 51 interface operator(/=) 52 module procedure class_ne 53 module procedure round_ne 54 end interface operator(/=) 55 56 ! See Fortran 2018, 17.10 & 17.11 57 interface ieee_class 58 module procedure ieee_class_a2 59 module procedure ieee_class_a3 60 module procedure ieee_class_a4 61 module procedure ieee_class_a8 62 module procedure ieee_class_a10 63 module procedure ieee_class_a16 64 end interface ieee_class 65 66 interface ieee_copy_sign 67 module procedure ieee_copy_sign_a2 68 module procedure ieee_copy_sign_a3 69 module procedure ieee_copy_sign_a4 70 module procedure ieee_copy_sign_a8 71 module procedure ieee_copy_sign_a10 72 module procedure ieee_copy_sign_a16 73 end interface ieee_copy_sign 74 75 ! TODO: more interfaces (_fma, &c.) 76 77 private :: classify 78 79 contains 80 81 elemental logical function class_eq(x,y) 82 type(ieee_class_type), intent(in) :: x, y 83 class_eq = x%which == y%which 84 end function class_eq 85 86 elemental logical function class_ne(x,y) 87 type(ieee_class_type), intent(in) :: x, y 88 class_ne = x%which /= y%which 89 end function class_ne 90 91 elemental logical function round_eq(x,y) 92 type(ieee_round_type), intent(in) :: x, y 93 round_eq = x%mode == y%mode 94 end function round_eq 95 96 elemental logical function round_ne(x,y) 97 type(ieee_round_type), intent(in) :: x, y 98 round_ne = x%mode /= y%mode 99 end function round_ne 100 101 elemental type(ieee_class_type) function classify( & 102 expo,maxExpo,negative,significandNZ,quietBit) 103 integer, intent(in) :: expo, maxExpo 104 logical, intent(in) :: negative, significandNZ, quietBit 105 if (expo == 0) then 106 if (significandNZ) then 107 if (negative) then 108 classify = ieee_negative_denormal 109 else 110 classify = ieee_positive_denormal 111 end if 112 else 113 if (negative) then 114 classify = ieee_negative_zero 115 else 116 classify = ieee_positive_zero 117 end if 118 end if 119 else if (expo == maxExpo) then 120 if (significandNZ) then 121 if (quietBit) then 122 classify = ieee_quiet_nan 123 else 124 classify = ieee_signaling_nan 125 end if 126 else 127 if (negative) then 128 classify = ieee_negative_inf 129 else 130 classify = ieee_positive_inf 131 end if 132 end if 133 else 134 if (negative) then 135 classify = ieee_negative_normal 136 else 137 classify = ieee_positive_normal 138 end if 139 end if 140 end function classify 141 142#define _CLASSIFY(RKIND,IKIND,TOTALBITS,PREC,IMPLICIT) \ 143 type(ieee_class_type) elemental function ieee_class_a##RKIND(x); \ 144 real(kind=RKIND), intent(in) :: x; \ 145 integer(kind=IKIND) :: raw; \ 146 integer, parameter :: significand = PREC - IMPLICIT; \ 147 integer, parameter :: exponentBits = TOTALBITS - 1 - significand; \ 148 integer, parameter :: maxExpo = shiftl(1, exponentBits) - 1; \ 149 integer :: exponent, sign; \ 150 logical :: negative, nzSignificand, quiet; \ 151 raw = transfer(x, raw); \ 152 exponent = ibits(raw, significand, exponentBits); \ 153 negative = btest(raw, TOTALBITS - 1); \ 154 nzSignificand = ibits(raw, 0, significand) /= 0; \ 155 quiet = btest(raw, significand - 1); \ 156 ieee_class_a##RKIND = classify(exponent, maxExpo, negative, nzSignificand, quiet); \ 157 end function ieee_class_a##RKIND 158 _CLASSIFY(2,2,16,11,1) 159 _CLASSIFY(3,2,16,8,1) 160 _CLASSIFY(4,4,32,24,1) 161 _CLASSIFY(8,8,64,53,1) 162 _CLASSIFY(10,16,80,64,0) 163 _CLASSIFY(16,16,128,112,1) 164#undef _CLASSIFY 165 166 ! TODO: This might need to be an actual Operation instead 167#define _COPYSIGN(RKIND,IKIND,BITS) \ 168 real(kind=RKIND) elemental function ieee_copy_sign_a##RKIND(x,y); \ 169 real(kind=RKIND), intent(in) :: x, y; \ 170 integer(kind=IKIND) :: xbits, ybits; \ 171 xbits = transfer(x, xbits); \ 172 ybits = transfer(y, ybits); \ 173 xbits = ior(ibclr(xbits, BITS-1), iand(ybits, shiftl(1_##IKIND, BITS-1))); \ 174 ieee_copy_sign_a##RKIND = transfer(xbits, x); \ 175 end function ieee_copy_sign_a##RKIND 176 _COPYSIGN(2,2,16) 177 _COPYSIGN(3,2,16) 178 _COPYSIGN(4,4,32) 179 _COPYSIGN(8,8,64) 180 _COPYSIGN(10,16,80) 181 _COPYSIGN(16,16,128) 182#undef _COPYSIGN 183 184end module ieee_arithmetic 185