xref: /aosp_15_r20/external/llvm-libc/src/math/generic/fmul.cpp (revision 71db0c75aadcf003ffe3238005f61d7618a3fead)
1 //===-- Implementation of fmul function -----------------------------------===//
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 #include "src/math/fmul.h"
9 #include "hdr/errno_macros.h"
10 #include "hdr/fenv_macros.h"
11 #include "src/__support/FPUtil/double_double.h"
12 #include "src/__support/FPUtil/generic/mul.h"
13 #include "src/__support/common.h"
14 #include "src/__support/macros/config.h"
15 
16 namespace LIBC_NAMESPACE_DECL {
17 
18 LLVM_LIBC_FUNCTION(float, fmul, (double x, double y)) {
19 
20   // Without FMA instructions, fputil::exact_mult is not
21   // correctly rounded for all rounding modes, so we fall
22   // back to the generic `fmul` implementation
23 
24 #ifndef LIBC_TARGET_CPU_HAS_FMA
25   return fputil::generic::mul<float>(x, y);
26 #else
27   fputil::DoubleDouble prod = fputil::exact_mult(x, y);
28   using DoubleBits = fputil::FPBits<double>;
29   using DoubleStorageType = typename DoubleBits::StorageType;
30   using FloatBits = fputil::FPBits<float>;
31   using FloatStorageType = typename FloatBits::StorageType;
32   DoubleBits x_bits(x);
33   DoubleBits y_bits(y);
34 
35   Sign result_sign = x_bits.sign() == y_bits.sign() ? Sign::POS : Sign::NEG;
36   double result = prod.hi;
37   DoubleBits hi_bits(prod.hi), lo_bits(prod.lo);
38   // Check for cases where we need to propagate the sticky bits:
39   constexpr uint64_t STICKY_MASK = 0xFFF'FFF; // Lower (52 - 23 - 1 = 28 bits)
40   uint64_t sticky_bits = (hi_bits.uintval() & STICKY_MASK);
41   if (LIBC_UNLIKELY(sticky_bits == 0)) {
42     // Might need to propagate sticky bits:
43     if (!(lo_bits.is_inf_or_nan() || lo_bits.is_zero())) {
44       // Now prod.lo is nonzero and finite, we need to propagate sticky bits.
45       if (lo_bits.sign() != hi_bits.sign())
46         result = DoubleBits(hi_bits.uintval() - 1).get_val();
47       else
48         result = DoubleBits(hi_bits.uintval() | 1).get_val();
49     }
50   }
51 
52   float result_f = static_cast<float>(result);
53   FloatBits rf_bits(result_f);
54   uint32_t rf_exp = rf_bits.get_biased_exponent();
55   if (LIBC_LIKELY(rf_exp > 0 && rf_exp < 2 * FloatBits::EXP_BIAS + 1)) {
56     return result_f;
57   }
58 
59   // Now result_f is either inf/nan/zero/denormal.
60   if (x_bits.is_nan() || y_bits.is_nan()) {
61     if (x_bits.is_signaling_nan() || y_bits.is_signaling_nan())
62       fputil::raise_except_if_required(FE_INVALID);
63 
64     if (x_bits.is_quiet_nan()) {
65       DoubleStorageType x_payload = x_bits.get_mantissa();
66       x_payload >>= DoubleBits::FRACTION_LEN - FloatBits::FRACTION_LEN;
67       return FloatBits::quiet_nan(x_bits.sign(),
68                                   static_cast<FloatStorageType>(x_payload))
69           .get_val();
70     }
71 
72     if (y_bits.is_quiet_nan()) {
73       DoubleStorageType y_payload = y_bits.get_mantissa();
74       y_payload >>= DoubleBits::FRACTION_LEN - FloatBits::FRACTION_LEN;
75       return FloatBits::quiet_nan(y_bits.sign(),
76                                   static_cast<FloatStorageType>(y_payload))
77           .get_val();
78     }
79 
80     return FloatBits::quiet_nan().get_val();
81   }
82 
83   if (x_bits.is_inf()) {
84     if (y_bits.is_zero()) {
85       fputil::set_errno_if_required(EDOM);
86       fputil::raise_except_if_required(FE_INVALID);
87 
88       return FloatBits::quiet_nan().get_val();
89     }
90 
91     return FloatBits::inf(result_sign).get_val();
92   }
93 
94   if (y_bits.is_inf()) {
95     if (x_bits.is_zero()) {
96       fputil::set_errno_if_required(EDOM);
97       fputil::raise_except_if_required(FE_INVALID);
98       return FloatBits::quiet_nan().get_val();
99     }
100 
101     return FloatBits::inf(result_sign).get_val();
102   }
103 
104   // Now either x or y is zero, and the other one is finite.
105   if (rf_bits.is_inf()) {
106     fputil::set_errno_if_required(ERANGE);
107     return FloatBits::inf(result_sign).get_val();
108   }
109 
110   if (x_bits.is_zero() || y_bits.is_zero())
111     return FloatBits::zero(result_sign).get_val();
112 
113   fputil::set_errno_if_required(ERANGE);
114   fputil::raise_except_if_required(FE_UNDERFLOW);
115   return result_f;
116 
117 #endif
118 }
119 } // namespace LIBC_NAMESPACE_DECL
120