73 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H11 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. We19// 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 use23// 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 static 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^102230 int64_t exp_hi = static_cast<int64_t>(hi + 1022) << FPBits::FRACTION_LEN;31 double mid_hi = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(mid));32 double lo_scaled =33 (lo != 0.0) ? cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(lo))34 : 0.0;35 36 double extra_factor = 0.0;37 uint64_t scale_down = 0x3FE0'0000'0000'0000; // 1022 in the exponent field.38 39 // Result is denormal if (mid_hi + lo_scale < 1.0).40 if ((1.0 - mid_hi) > lo_scaled) {41 // Extra rounding step is needed, which adds more rounding errors.42 err += 0x1.0p-52;43 extra_factor = 1.0;44 scale_down = 0x3FF0'0000'0000'0000; // 1023 in the exponent field.45 }46 47 // By adding 1.0, the results will have similar rounding points as denormal48 // outputs.49 if constexpr (SKIP_ZIV_TEST) {50 double r = extra_factor + (mid_hi + lo_scaled);51 return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(r) - scale_down);52 } else {53 double err_scaled =54 cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(err));55 56 double lo_u = lo_scaled + err_scaled;57 double lo_l = lo_scaled - err_scaled;58 59 double upper = extra_factor + (mid_hi + lo_u);60 double lower = extra_factor + (mid_hi + lo_l);61 62 if (LIBC_LIKELY(upper == lower)) {63 return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(upper) - scale_down);64 }65 66 return cpp::nullopt;67 }68}69 70} // namespace LIBC_NAMESPACE_DECL71 72#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H73