1 //===-- A class to store a normalized floating point number -----*- C++ -*-===//
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 #ifndef LLVM_LIBC_SRC_SUPPORT_FPUTIL_NORMAL_FLOAT_H
10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_NORMAL_FLOAT_H
11 
12 #include "FPBits.h"
13 
14 #include "src/__support/CPP/TypeTraits.h"
15 
16 #include <stdint.h>
17 
18 namespace __llvm_libc {
19 namespace fputil {
20 
21 // A class which stores the normalized form of a floating point value.
22 // The special IEEE-754 bits patterns of Zero, infinity and NaNs are
23 // are not handled by this class.
24 //
25 // A normalized floating point number is of this form:
26 //    (-1)*sign * 2^exponent * <mantissa>
27 // where <mantissa> is of the form 1.<...>.
28 template <typename T> struct NormalFloat {
29   static_assert(
30       cpp::IsFloatingPointType<T>::Value,
31       "NormalFloat template parameter has to be a floating point type.");
32 
33   using UIntType = typename FPBits<T>::UIntType;
34   static constexpr UIntType ONE = (UIntType(1) << MantissaWidth<T>::VALUE);
35 
36   // Unbiased exponent value.
37   int32_t exponent;
38 
39   UIntType mantissa;
40   // We want |UIntType| to have atleast one bit more than the actual mantissa
41   // bit width to accommodate the implicit 1 value.
42   static_assert(sizeof(UIntType) * 8 >= MantissaWidth<T>::VALUE + 1,
43                 "Bad type for mantissa in NormalFloat.");
44 
45   bool sign;
46 
47   NormalFloat(int32_t e, UIntType m, bool s)
48       : exponent(e), mantissa(m), sign(s) {
49     if (mantissa >= ONE)
50       return;
51 
52     unsigned normalization_shift = evaluate_normalization_shift(mantissa);
53     mantissa = mantissa << normalization_shift;
54     exponent -= normalization_shift;
55   }
56 
57   explicit NormalFloat(T x) { init_from_bits(FPBits<T>(x)); }
58 
59   explicit NormalFloat(FPBits<T> bits) { init_from_bits(bits); }
60 
61   // Compares this normalized number with another normalized number.
62   // Returns -1 is this number is less than |other|, 0 if this number is equal
63   // to |other|, and 1 if this number is greater than |other|.
64   int cmp(const NormalFloat<T> &other) const {
65     if (sign != other.sign)
66       return sign ? -1 : 1;
67 
68     if (exponent > other.exponent) {
69       return sign ? -1 : 1;
70     } else if (exponent == other.exponent) {
71       if (mantissa > other.mantissa)
72         return sign ? -1 : 1;
73       else if (mantissa == other.mantissa)
74         return 0;
75       else
76         return sign ? 1 : -1;
77     } else {
78       return sign ? 1 : -1;
79     }
80   }
81 
82   // Returns a new normalized floating point number which is equal in value
83   // to this number multiplied by 2^e. That is:
84   //     new = this *  2^e
85   NormalFloat<T> mul2(int e) const {
86     NormalFloat<T> result = *this;
87     result.exponent += e;
88     return result;
89   }
90 
91   operator T() const {
92     int biased_exponent = exponent + FPBits<T>::EXPONENT_BIAS;
93     // Max exponent is of the form 0xFF...E. That is why -2 and not -1.
94     constexpr int MAX_EXPONENT_VALUE = (1 << ExponentWidth<T>::VALUE) - 2;
95     if (biased_exponent > MAX_EXPONENT_VALUE) {
96       return sign ? T(FPBits<T>::neg_inf()) : T(FPBits<T>::inf());
97     }
98 
99     FPBits<T> result(T(0.0));
100     result.set_sign(sign);
101 
102     constexpr int SUBNORMAL_EXPONENT = -FPBits<T>::EXPONENT_BIAS + 1;
103     if (exponent < SUBNORMAL_EXPONENT) {
104       unsigned shift = SUBNORMAL_EXPONENT - exponent;
105       // Since exponent > subnormalExponent, shift is strictly greater than
106       // zero.
107       if (shift <= MantissaWidth<T>::VALUE + 1) {
108         // Generate a subnormal number. Might lead to loss of precision.
109         // We round to nearest and round halfway cases to even.
110         const UIntType shift_out_mask = (UIntType(1) << shift) - 1;
111         const UIntType shift_out_value = mantissa & shift_out_mask;
112         const UIntType halfway_value = UIntType(1) << (shift - 1);
113         result.set_unbiased_exponent(0);
114         result.set_mantissa(mantissa >> shift);
115         UIntType new_mantissa = result.get_mantissa();
116         if (shift_out_value > halfway_value) {
117           new_mantissa += 1;
118         } else if (shift_out_value == halfway_value) {
119           // Round to even.
120           if (result.get_mantissa() & 0x1)
121             new_mantissa += 1;
122         }
123         result.set_mantissa(new_mantissa);
124         // Adding 1 to mantissa can lead to overflow. This can only happen if
125         // mantissa was all ones (0b111..11). For such a case, we will carry
126         // the overflow into the exponent.
127         if (new_mantissa == ONE)
128           result.set_unbiased_exponent(1);
129         return T(result);
130       } else {
131         return T(result);
132       }
133     }
134 
135     result.set_unbiased_exponent(exponent + FPBits<T>::EXPONENT_BIAS);
136     result.set_mantissa(mantissa);
137     return T(result);
138   }
139 
140 private:
141   void init_from_bits(FPBits<T> bits) {
142     sign = bits.get_sign();
143 
144     if (bits.is_inf_or_nan() || bits.is_zero()) {
145       // Ignore special bit patterns. Implementations deal with them separately
146       // anyway so this should not be a problem.
147       exponent = 0;
148       mantissa = 0;
149       return;
150     }
151 
152     // Normalize subnormal numbers.
153     if (bits.get_unbiased_exponent() == 0) {
154       unsigned shift = evaluate_normalization_shift(bits.get_mantissa());
155       mantissa = UIntType(bits.get_mantissa()) << shift;
156       exponent = 1 - FPBits<T>::EXPONENT_BIAS - shift;
157     } else {
158       exponent = bits.get_unbiased_exponent() - FPBits<T>::EXPONENT_BIAS;
159       mantissa = ONE | bits.get_mantissa();
160     }
161   }
162 
163   unsigned evaluate_normalization_shift(UIntType m) {
164     unsigned shift = 0;
165     for (; (ONE & m) == 0 && (shift < MantissaWidth<T>::VALUE);
166          m <<= 1, ++shift)
167       ;
168     return shift;
169   }
170 };
171 
172 #ifdef SPECIAL_X86_LONG_DOUBLE
173 template <>
174 inline void NormalFloat<long double>::init_from_bits(FPBits<long double> bits) {
175   sign = bits.get_sign();
176 
177   if (bits.is_inf_or_nan() || bits.is_zero()) {
178     // Ignore special bit patterns. Implementations deal with them separately
179     // anyway so this should not be a problem.
180     exponent = 0;
181     mantissa = 0;
182     return;
183   }
184 
185   if (bits.get_unbiased_exponent() == 0) {
186     if (bits.get_implicit_bit() == 0) {
187       // Since we ignore zero value, the mantissa in this case is non-zero.
188       int normalization_shift =
189           evaluate_normalization_shift(bits.get_mantissa());
190       exponent = -16382 - normalization_shift;
191       mantissa = (bits.get_mantissa() << normalization_shift);
192     } else {
193       exponent = -16382;
194       mantissa = ONE | bits.get_mantissa();
195     }
196   } else {
197     if (bits.get_implicit_bit() == 0) {
198       // Invalid number so just store 0 similar to a NaN.
199       exponent = 0;
200       mantissa = 0;
201     } else {
202       exponent = bits.get_unbiased_exponent() - 16383;
203       mantissa = ONE | bits.get_mantissa();
204     }
205   }
206 }
207 
208 template <> inline NormalFloat<long double>::operator long double() const {
209   int biased_exponent = exponent + FPBits<long double>::EXPONENT_BIAS;
210   // Max exponent is of the form 0xFF...E. That is why -2 and not -1.
211   constexpr int MAX_EXPONENT_VALUE =
212       (1 << ExponentWidth<long double>::VALUE) - 2;
213   if (biased_exponent > MAX_EXPONENT_VALUE) {
214     return sign ? FPBits<long double>::neg_inf() : FPBits<long double>::inf();
215   }
216 
217   FPBits<long double> result(0.0l);
218   result.set_sign(sign);
219 
220   constexpr int SUBNORMAL_EXPONENT = -FPBits<long double>::EXPONENT_BIAS + 1;
221   if (exponent < SUBNORMAL_EXPONENT) {
222     unsigned shift = SUBNORMAL_EXPONENT - exponent;
223     if (shift <= MantissaWidth<long double>::VALUE + 1) {
224       // Generate a subnormal number. Might lead to loss of precision.
225       // We round to nearest and round halfway cases to even.
226       const UIntType shift_out_mask = (UIntType(1) << shift) - 1;
227       const UIntType shift_out_value = mantissa & shift_out_mask;
228       const UIntType halfway_value = UIntType(1) << (shift - 1);
229       result.set_unbiased_exponent(0);
230       result.set_mantissa(mantissa >> shift);
231       UIntType new_mantissa = result.get_mantissa();
232       if (shift_out_value > halfway_value) {
233         new_mantissa += 1;
234       } else if (shift_out_value == halfway_value) {
235         // Round to even.
236         if (result.get_mantissa() & 0x1)
237           new_mantissa += 1;
238       }
239       result.set_mantissa(new_mantissa);
240       // Adding 1 to mantissa can lead to overflow. This can only happen if
241       // mantissa was all ones (0b111..11). For such a case, we will carry
242       // the overflow into the exponent and set the implicit bit to 1.
243       if (new_mantissa == ONE) {
244         result.set_unbiased_exponent(1);
245         result.set_implicit_bit(1);
246       } else {
247         result.set_implicit_bit(0);
248       }
249       return static_cast<long double>(result);
250     } else {
251       return static_cast<long double>(result);
252     }
253   }
254 
255   result.set_unbiased_exponent(biased_exponent);
256   result.set_mantissa(mantissa);
257   result.set_implicit_bit(1);
258   return static_cast<long double>(result);
259 }
260 #endif // SPECIAL_X86_LONG_DOUBLE
261 
262 } // namespace fputil
263 } // namespace __llvm_libc
264 
265 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_NORMAL_FLOAT_H
266