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