1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
|
//===-- 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_DOUBLE
return fputil::generic::mul<float>(x, y);
#else
fputil::DoubleDouble prod = fputil::exact_mult(x, y);
using DoubleBits = fputil::FPBits<double>;
using DoubleStorageType = typename DoubleBits::StorageType;
using FloatBits = fputil::FPBits<float>;
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<float>(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<FloatStorageType>(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<FloatStorageType>(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
|
