120 lines · cpp
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-exception6//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 not21 // correctly rounded for all rounding modes, so we fall22 // back to the generic `fmul` implementation23 24#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE25 return fputil::generic::mul<float>(x, y);26#else27 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 else48 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#endif118}119} // namespace LIBC_NAMESPACE_DECL120