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
16namespace LIBC_NAMESPACE_DECL {
17
18LLVM_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_DOUBLE
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

source code of libc/src/math/generic/fmul.cpp