1//===-- Common utils for exp function ---------------------------*- 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_EXP_UTILS_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
11
12#include "src/__support/CPP/bit.h"
13#include "src/__support/CPP/optional.h"
14#include "src/__support/FPUtil/FPBits.h"
15
16namespace LIBC_NAMESPACE_DECL {
17
18// Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We
19// assume further that 1 <= mid < 2, mid + lo < 2, and |lo| << mid.
20// Notice that, if 0 < x < 2^-1022,
21// double(2^-1022 + x) - 2^-1022 = double(x).
22// So if we scale x up by 2^1022, we can use
23// double(1.0 + 2^1022 * x) - 1.0 to test how x is rounded in denormal range.
24template <bool SKIP_ZIV_TEST = false>
25LIBC_INLINE constexpr cpp::optional<double>
26ziv_test_denorm(int hi, double mid, double lo, double err) {
27 using FPBits = typename fputil::FPBits<double>;
28
29 // Scaling factor = 1/(min normal number) = 2^1022
30 int64_t exp_hi =
31 static_cast<int64_t>(hi + 1022) * (1LL << FPBits::FRACTION_LEN);
32 double mid_hi = cpp::bit_cast<double>(from: exp_hi + cpp::bit_cast<int64_t>(from: mid));
33 double lo_scaled =
34 (lo != 0.0) ? cpp::bit_cast<double>(from: exp_hi + cpp::bit_cast<int64_t>(from: lo))
35 : 0.0;
36
37 double extra_factor = 0.0;
38 uint64_t scale_down = 0x3FE0'0000'0000'0000; // 1022 in the exponent field.
39
40 // Result is denormal if (mid_hi + lo_scale < 1.0).
41 if ((1.0 - mid_hi) > lo_scaled) {
42 // Extra rounding step is needed, which adds more rounding errors.
43 err += 0x1.0p-52;
44 extra_factor = 1.0;
45 scale_down = 0x3FF0'0000'0000'0000; // 1023 in the exponent field.
46 }
47
48 // By adding 1.0, the results will have similar rounding points as denormal
49 // outputs.
50 if constexpr (SKIP_ZIV_TEST) {
51 double r = extra_factor + (mid_hi + lo_scaled);
52 return cpp::bit_cast<double>(from: cpp::bit_cast<uint64_t>(from: r) - scale_down);
53 } else {
54 double err_scaled =
55 cpp::bit_cast<double>(from: exp_hi + cpp::bit_cast<int64_t>(from: err));
56
57 double lo_u = lo_scaled + err_scaled;
58 double lo_l = lo_scaled - err_scaled;
59
60 double upper = extra_factor + (mid_hi + lo_u);
61 double lower = extra_factor + (mid_hi + lo_l);
62
63 if (LIBC_LIKELY(upper == lower)) {
64 return cpp::bit_cast<double>(from: cpp::bit_cast<uint64_t>(from: upper) - scale_down);
65 }
66
67 return cpp::nullopt;
68 }
69}
70
71} // namespace LIBC_NAMESPACE_DECL
72
73#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
74