//===-- Implementation of fmul function -----------------------------------===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// #include "src/math/fmul.h" #include "hdr/errno_macros.h" #include "hdr/fenv_macros.h" #include "src/__support/FPUtil/double_double.h" #include "src/__support/FPUtil/generic/mul.h" #include "src/__support/common.h" #include "src/__support/macros/config.h" namespace LIBC_NAMESPACE_DECL { LLVM_LIBC_FUNCTION(float, fmul, (double x, double y)) { // Without FMA instructions, fputil::exact_mult is not // correctly rounded for all rounding modes, so we fall // back to the generic `fmul` implementation #ifndef LIBC_TARGET_CPU_HAS_FMA return fputil::generic::mul(x, y); #else fputil::DoubleDouble prod = fputil::exact_mult(x, y); using DoubleBits = fputil::FPBits; using DoubleStorageType = typename DoubleBits::StorageType; using FloatBits = fputil::FPBits; using FloatStorageType = typename FloatBits::StorageType; DoubleBits x_bits(x); DoubleBits y_bits(y); Sign result_sign = x_bits.sign() == y_bits.sign() ? Sign::POS : Sign::NEG; double result = prod.hi; DoubleBits hi_bits(prod.hi), lo_bits(prod.lo); // Check for cases where we need to propagate the sticky bits: constexpr uint64_t STICKY_MASK = 0xFFF'FFF; // Lower (52 - 23 - 1 = 28 bits) uint64_t sticky_bits = (hi_bits.uintval() & STICKY_MASK); if (LIBC_UNLIKELY(sticky_bits == 0)) { // Might need to propagate sticky bits: if (!(lo_bits.is_inf_or_nan() || lo_bits.is_zero())) { // Now prod.lo is nonzero and finite, we need to propagate sticky bits. if (lo_bits.sign() != hi_bits.sign()) result = DoubleBits(hi_bits.uintval() - 1).get_val(); else result = DoubleBits(hi_bits.uintval() | 1).get_val(); } } float result_f = static_cast(result); FloatBits rf_bits(result_f); uint32_t rf_exp = rf_bits.get_biased_exponent(); if (LIBC_LIKELY(rf_exp > 0 && rf_exp < 2 * FloatBits::EXP_BIAS + 1)) { return result_f; } // Now result_f is either inf/nan/zero/denormal. if (x_bits.is_nan() || y_bits.is_nan()) { if (x_bits.is_signaling_nan() || y_bits.is_signaling_nan()) fputil::raise_except_if_required(FE_INVALID); if (x_bits.is_quiet_nan()) { DoubleStorageType x_payload = x_bits.get_mantissa(); x_payload >>= DoubleBits::FRACTION_LEN - FloatBits::FRACTION_LEN; return FloatBits::quiet_nan(x_bits.sign(), static_cast(x_payload)) .get_val(); } if (y_bits.is_quiet_nan()) { DoubleStorageType y_payload = y_bits.get_mantissa(); y_payload >>= DoubleBits::FRACTION_LEN - FloatBits::FRACTION_LEN; return FloatBits::quiet_nan(y_bits.sign(), static_cast(y_payload)) .get_val(); } return FloatBits::quiet_nan().get_val(); } if (x_bits.is_inf()) { if (y_bits.is_zero()) { fputil::set_errno_if_required(EDOM); fputil::raise_except_if_required(FE_INVALID); return FloatBits::quiet_nan().get_val(); } return FloatBits::inf(result_sign).get_val(); } if (y_bits.is_inf()) { if (x_bits.is_zero()) { fputil::set_errno_if_required(EDOM); fputil::raise_except_if_required(FE_INVALID); return FloatBits::quiet_nan().get_val(); } return FloatBits::inf(result_sign).get_val(); } // Now either x or y is zero, and the other one is finite. if (rf_bits.is_inf()) { fputil::set_errno_if_required(ERANGE); return FloatBits::inf(result_sign).get_val(); } if (x_bits.is_zero() || y_bits.is_zero()) return FloatBits::zero(result_sign).get_val(); fputil::set_errno_if_required(ERANGE); fputil::raise_except_if_required(FE_UNDERFLOW); return result_f; #endif } } // namespace LIBC_NAMESPACE_DECL