1 //===-- Common header for FMA implementations -------------------*- 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_GENERIC_FMA_H
10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMA_H
11 
12 #include "src/__support/CPP/Bit.h"
13 #include "src/__support/CPP/TypeTraits.h"
14 #include "src/__support/CPP/UInt128.h"
15 #include "src/__support/FPUtil/FEnvImpl.h"
16 #include "src/__support/FPUtil/FPBits.h"
17 #include "src/__support/FPUtil/FloatProperties.h"
18 #include "src/__support/FPUtil/builtin_wrappers.h"
19 #include "src/__support/common.h"
20 
21 namespace __llvm_libc {
22 namespace fputil {
23 namespace generic {
24 
25 template <typename T> static inline T fma(T x, T y, T z);
26 
27 // TODO(lntue): Implement fmaf that is correctly rounded to all rounding modes.
28 // The implementation below only is only correct for the default rounding mode,
29 // round-to-nearest tie-to-even.
30 template <> inline float fma<float>(float x, float y, float z) {
31   // Product is exact.
32   double prod = static_cast<double>(x) * static_cast<double>(y);
33   double z_d = static_cast<double>(z);
34   double sum = prod + z_d;
35   fputil::FPBits<double> bit_prod(prod), bitz(z_d), bit_sum(sum);
36 
37   if (!(bit_sum.is_inf_or_nan() || bit_sum.is_zero())) {
38     // Since the sum is computed in double precision, rounding might happen
39     // (for instance, when bitz.exponent > bit_prod.exponent + 5, or
40     // bit_prod.exponent > bitz.exponent + 40).  In that case, when we round
41     // the sum back to float, double rounding error might occur.
42     // A concrete example of this phenomenon is as follows:
43     //   x = y = 1 + 2^(-12), z = 2^(-53)
44     // The exact value of x*y + z is 1 + 2^(-11) + 2^(-24) + 2^(-53)
45     // So when rounding to float, fmaf(x, y, z) = 1 + 2^(-11) + 2^(-23)
46     // On the other hand, with the default rounding mode,
47     //   double(x*y + z) = 1 + 2^(-11) + 2^(-24)
48     // and casting again to float gives us:
49     //   float(double(x*y + z)) = 1 + 2^(-11).
50     //
51     // In order to correct this possible double rounding error, first we use
52     // Dekker's 2Sum algorithm to find t such that sum - t = prod + z exactly,
53     // assuming the (default) rounding mode is round-to-the-nearest,
54     // tie-to-even.  Moreover, t satisfies the condition that t < eps(sum),
55     // i.e., t.exponent < sum.exponent - 52. So if t is not 0, meaning rounding
56     // occurs when computing the sum, we just need to use t to adjust (any) last
57     // bit of sum, so that the sticky bits used when rounding sum to float are
58     // correct (when it matters).
59     fputil::FPBits<double> t(
60         (bit_prod.get_unbiased_exponent() >= bitz.get_unbiased_exponent())
61             ? ((double(bit_sum) - double(bit_prod)) - double(bitz))
62             : ((double(bit_sum) - double(bitz)) - double(bit_prod)));
63 
64     // Update sticky bits if t != 0.0 and the least (52 - 23 - 1 = 28) bits are
65     // zero.
66     if (!t.is_zero() && ((bit_sum.get_mantissa() & 0xfff'ffffULL) == 0)) {
67       if (bit_sum.get_sign() != t.get_sign()) {
68         bit_sum.set_mantissa(bit_sum.get_mantissa() + 1);
69       } else if (bit_sum.get_mantissa()) {
70         bit_sum.set_mantissa(bit_sum.get_mantissa() - 1);
71       }
72     }
73   }
74 
75   return static_cast<float>(static_cast<double>(bit_sum));
76 }
77 
78 namespace internal {
79 
80 // Extract the sticky bits and shift the `mantissa` to the right by
81 // `shift_length`.
82 static inline bool shift_mantissa(int shift_length, UInt128 &mant) {
83   if (shift_length >= 128) {
84     mant = 0;
85     return true; // prod_mant is non-zero.
86   }
87   UInt128 mask = (UInt128(1) << shift_length) - 1;
88   bool sticky_bits = (mant & mask) != 0;
89   mant >>= shift_length;
90   return sticky_bits;
91 }
92 
93 } // namespace internal
94 
95 template <> inline double fma<double>(double x, double y, double z) {
96   using FPBits = fputil::FPBits<double>;
97   using FloatProp = fputil::FloatProperties<double>;
98 
99   if (unlikely(x == 0 || y == 0 || z == 0)) {
100     return x * y + z;
101   }
102 
103   int x_exp = 0;
104   int y_exp = 0;
105   int z_exp = 0;
106 
107   // Normalize denormal inputs.
108   if (unlikely(FPBits(x).get_unbiased_exponent() == 0)) {
109     x_exp -= 52;
110     x *= 0x1.0p+52;
111   }
112   if (unlikely(FPBits(y).get_unbiased_exponent() == 0)) {
113     y_exp -= 52;
114     y *= 0x1.0p+52;
115   }
116   if (unlikely(FPBits(z).get_unbiased_exponent() == 0)) {
117     z_exp -= 52;
118     z *= 0x1.0p+52;
119   }
120 
121   FPBits x_bits(x), y_bits(y), z_bits(z);
122   bool x_sign = x_bits.get_sign();
123   bool y_sign = y_bits.get_sign();
124   bool z_sign = z_bits.get_sign();
125   bool prod_sign = x_sign != y_sign;
126   x_exp += x_bits.get_unbiased_exponent();
127   y_exp += y_bits.get_unbiased_exponent();
128   z_exp += z_bits.get_unbiased_exponent();
129 
130   if (unlikely(x_exp == FPBits::MAX_EXPONENT || y_exp == FPBits::MAX_EXPONENT ||
131                z_exp == FPBits::MAX_EXPONENT))
132     return x * y + z;
133 
134   // Extract mantissa and append hidden leading bits.
135   UInt128 x_mant = x_bits.get_mantissa() | FPBits::MIN_NORMAL;
136   UInt128 y_mant = y_bits.get_mantissa() | FPBits::MIN_NORMAL;
137   UInt128 z_mant = z_bits.get_mantissa() | FPBits::MIN_NORMAL;
138 
139   // If the exponent of the product x*y > the exponent of z, then no extra
140   // precision beside the entire product x*y is needed.  On the other hand, when
141   // the exponent of z >= the exponent of the product x*y, the worst-case that
142   // we need extra precision is when there is cancellation and the most
143   // significant bit of the product is aligned exactly with the second most
144   // significant bit of z:
145   //      z :    10aa...a
146   // - prod :     1bb...bb....b
147   // In that case, in order to store the exact result, we need at least
148   //   (Length of prod) - (MantissaLength of z) = 2*(52 + 1) - 52 = 54.
149   // Overall, before aligning the mantissas and exponents, we can simply left-
150   // shift the mantissa of z by at least 54, and left-shift the product of x*y
151   // by (that amount - 52).  After that, it is enough to align the least
152   // significant bit, given that we keep track of the round and sticky bits
153   // after the least significant bit.
154   // We pick shifting z_mant by 64 bits so that technically we can simply use
155   // the original mantissa as high part when constructing 128-bit z_mant. So the
156   // mantissa of prod will be left-shifted by 64 - 54 = 10 initially.
157 
158   UInt128 prod_mant = x_mant * y_mant << 10;
159   int prod_lsb_exp =
160       x_exp + y_exp -
161       (FPBits::EXPONENT_BIAS + 2 * MantissaWidth<double>::VALUE + 10);
162 
163   z_mant <<= 64;
164   int z_lsb_exp = z_exp - (MantissaWidth<double>::VALUE + 64);
165   bool round_bit = false;
166   bool sticky_bits = false;
167   bool z_shifted = false;
168 
169   // Align exponents.
170   if (prod_lsb_exp < z_lsb_exp) {
171     sticky_bits = internal::shift_mantissa(z_lsb_exp - prod_lsb_exp, prod_mant);
172     prod_lsb_exp = z_lsb_exp;
173   } else if (z_lsb_exp < prod_lsb_exp) {
174     z_shifted = true;
175     sticky_bits = internal::shift_mantissa(prod_lsb_exp - z_lsb_exp, z_mant);
176   }
177 
178   // Perform the addition:
179   //   (-1)^prod_sign * prod_mant + (-1)^z_sign * z_mant.
180   // The final result will be stored in prod_sign and prod_mant.
181   if (prod_sign == z_sign) {
182     // Effectively an addition.
183     prod_mant += z_mant;
184   } else {
185     // Subtraction cases.
186     if (prod_mant >= z_mant) {
187       if (z_shifted && sticky_bits) {
188         // Add 1 more to the subtrahend so that the sticky bits remain
189         // positive. This would simplify the rounding logic.
190         ++z_mant;
191       }
192       prod_mant -= z_mant;
193     } else {
194       if (!z_shifted && sticky_bits) {
195         // Add 1 more to the subtrahend so that the sticky bits remain
196         // positive. This would simplify the rounding logic.
197         ++prod_mant;
198       }
199       prod_mant = z_mant - prod_mant;
200       prod_sign = z_sign;
201     }
202   }
203 
204   uint64_t result = 0;
205   int r_exp = 0; // Unbiased exponent of the result
206 
207   // Normalize the result.
208   if (prod_mant != 0) {
209     uint64_t prod_hi = static_cast<uint64_t>(prod_mant >> 64);
210     int lead_zeros = prod_hi
211                          ? unsafe_clz(prod_hi)
212                          : 64 + unsafe_clz(static_cast<uint64_t>(prod_mant));
213     // Move the leading 1 to the most significant bit.
214     prod_mant <<= lead_zeros;
215     // The lower 64 bits are always sticky bits after moving the leading 1 to
216     // the most significant bit.
217     sticky_bits |= (static_cast<uint64_t>(prod_mant) != 0);
218     result = static_cast<uint64_t>(prod_mant >> 64);
219     // Change prod_lsb_exp the be the exponent of the least significant bit of
220     // the result.
221     prod_lsb_exp += 64 - lead_zeros;
222     r_exp = prod_lsb_exp + 63;
223 
224     if (r_exp > 0) {
225       // The result is normal.  We will shift the mantissa to the right by
226       // 63 - 52 = 11 bits (from the locations of the most significant bit).
227       // Then the rounding bit will correspond the the 11th bit, and the lowest
228       // 10 bits are merged into sticky bits.
229       round_bit = (result & 0x0400ULL) != 0;
230       sticky_bits |= (result & 0x03ffULL) != 0;
231       result >>= 11;
232     } else {
233       if (r_exp < -52) {
234         // The result is smaller than 1/2 of the smallest denormal number.
235         sticky_bits = true; // since the result is non-zero.
236         result = 0;
237       } else {
238         // The result is denormal.
239         uint64_t mask = 1ULL << (11 - r_exp);
240         round_bit = (result & mask) != 0;
241         sticky_bits |= (result & (mask - 1)) != 0;
242         if (r_exp > -52)
243           result >>= 12 - r_exp;
244         else
245           result = 0;
246       }
247 
248       r_exp = 0;
249     }
250   } else {
251     // Return +0.0 when there is exact cancellation, i.e., x*y == -z exactly.
252     prod_sign = false;
253   }
254 
255   // Finalize the result.
256   int round_mode = fputil::get_round();
257   if (unlikely(r_exp >= FPBits::MAX_EXPONENT)) {
258     if ((round_mode == FE_TOWARDZERO) ||
259         (round_mode == FE_UPWARD && prod_sign) ||
260         (round_mode == FE_DOWNWARD && !prod_sign)) {
261       result = FPBits::MAX_NORMAL;
262       return prod_sign ? -bit_cast<double>(result) : bit_cast<double>(result);
263     }
264     return prod_sign ? static_cast<double>(FPBits::neg_inf())
265                      : static_cast<double>(FPBits::inf());
266   }
267 
268   // Remove hidden bit and append the exponent field and sign bit.
269   result = (result & FloatProp::MANTISSA_MASK) |
270            (static_cast<uint64_t>(r_exp) << FloatProp::MANTISSA_WIDTH);
271   if (prod_sign) {
272     result |= FloatProp::SIGN_MASK;
273   }
274 
275   // Rounding.
276   if (round_mode == FE_TONEAREST) {
277     if (round_bit && (sticky_bits || ((result & 1) != 0)))
278       ++result;
279   } else if ((round_mode == FE_UPWARD && !prod_sign) ||
280              (round_mode == FE_DOWNWARD && prod_sign)) {
281     if (round_bit || sticky_bits)
282       ++result;
283   }
284 
285   return bit_cast<double>(result);
286 }
287 
288 } // namespace generic
289 } // namespace fputil
290 } // namespace __llvm_libc
291 
292 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMA_H
293