|
| 1 | +//===-- Implementation header for asinhf16 ----------------------*- 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_ASINHF16_H |
| 10 | +#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF16_H |
| 11 | + |
| 12 | +#include "include/llvm-libc-macros/float16-macros.h" |
| 13 | + |
| 14 | +#ifdef LIBC_TYPES_HAS_FLOAT16 |
| 15 | + |
| 16 | +#include "acoshf_utils.h" |
| 17 | +#include "src/__support/FPUtil/FEnvImpl.h" |
| 18 | +#include "src/__support/FPUtil/FPBits.h" |
| 19 | +#include "src/__support/FPUtil/PolyEval.h" |
| 20 | +#include "src/__support/FPUtil/cast.h" |
| 21 | +#include "src/__support/FPUtil/except_value_utils.h" |
| 22 | +#include "src/__support/FPUtil/multiply_add.h" |
| 23 | +#include "src/__support/FPUtil/rounding_mode.h" |
| 24 | +#include "src/__support/FPUtil/sqrt.h" |
| 25 | +#include "src/__support/macros/config.h" |
| 26 | +#include "src/__support/macros/optimization.h" |
| 27 | + |
| 28 | +namespace LIBC_NAMESPACE_DECL { |
| 29 | + |
| 30 | +namespace math { |
| 31 | + |
| 32 | +LIBC_INLINE static constexpr float16 asinhf16(float16 x) { |
| 33 | + |
| 34 | +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS |
| 35 | + constexpr size_t N_EXCEPTS = 8; |
| 36 | + |
| 37 | + constexpr fputil::ExceptValues<float16, N_EXCEPTS> ASINHF16_EXCEPTS{{ |
| 38 | + // (input, RZ output, RU offset, RD offset, RN offset) |
| 39 | + |
| 40 | + // x = 0x1.da4p-2, asinhf16(x) = 0x1.ca8p-2 (RZ) |
| 41 | + {0x3769, 0x372a, 1, 0, 1}, |
| 42 | + // x = 0x1.d6cp-1, asinhf16(x) = 0x1.a58p-1 (RZ) |
| 43 | + {0x3b5b, 0x3a96, 1, 0, 0}, |
| 44 | + // x = 0x1.c7cp+3, asinhf16(x) = 0x1.accp+1 (RZ) |
| 45 | + {0x4b1f, 0x42b3, 1, 0, 0}, |
| 46 | + // x = 0x1.26cp+4, asinhf16(x) = 0x1.cd8p+1 (RZ) |
| 47 | + {0x4c9b, 0x4336, 1, 0, 1}, |
| 48 | + // x = -0x1.da4p-2, asinhf16(x) = -0x1.ca8p-2 (RZ) |
| 49 | + {0xb769, 0xb72a, 0, 1, 1}, |
| 50 | + // x = -0x1.d6cp-1, asinhf16(x) = -0x1.a58p-1 (RZ) |
| 51 | + {0xbb5b, 0xba96, 0, 1, 0}, |
| 52 | + // x = -0x1.c7cp+3, asinhf16(x) = -0x1.accp+1 (RZ) |
| 53 | + {0xcb1f, 0xc2b3, 0, 1, 0}, |
| 54 | + // x = -0x1.26cp+4, asinhf16(x) = -0x1.cd8p+1 (RZ) |
| 55 | + {0xcc9b, 0xc336, 0, 1, 1}, |
| 56 | + }}; |
| 57 | +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS |
| 58 | + |
| 59 | + using namespace acoshf_internal; |
| 60 | + using FPBits = fputil::FPBits<float16>; |
| 61 | + FPBits xbits(x); |
| 62 | + |
| 63 | + uint16_t x_u = xbits.uintval(); |
| 64 | + uint16_t x_abs = x_u & 0x7fff; |
| 65 | + |
| 66 | + if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { |
| 67 | + if (xbits.is_signaling_nan()) { |
| 68 | + fputil::raise_except_if_required(FE_INVALID); |
| 69 | + return FPBits::quiet_nan().get_val(); |
| 70 | + } |
| 71 | + |
| 72 | + return x; |
| 73 | + } |
| 74 | + |
| 75 | +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS |
| 76 | + // Handle exceptional values |
| 77 | + if (auto r = ASINHF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value())) |
| 78 | + return r.value(); |
| 79 | +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS |
| 80 | + |
| 81 | + float xf = x; |
| 82 | + const float SIGN[2] = {1.0f, -1.0f}; |
| 83 | + float x_sign = SIGN[x_u >> 15]; |
| 84 | + |
| 85 | + // |x| <= 0.25 |
| 86 | + if (LIBC_UNLIKELY(x_abs <= 0x3400)) { |
| 87 | + // when |x| < 0x1.718p-5, asinhf16(x) = x. Adjust by 1 ULP for certain |
| 88 | + // rounding types. |
| 89 | + if (LIBC_UNLIKELY(x_abs < 0x29c6)) { |
| 90 | + int rounding = fputil::quick_get_round(); |
| 91 | + if ((rounding == FE_UPWARD || rounding == FE_TOWARDZERO) && xf < 0) |
| 92 | + return fputil::cast<float16>(xf + 0x1p-24f); |
| 93 | + if ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) && xf > 0) |
| 94 | + return fputil::cast<float16>(xf - 0x1p-24f); |
| 95 | + return fputil::cast<float16>(xf); |
| 96 | + } |
| 97 | + |
| 98 | + float x_sq = xf * xf; |
| 99 | + // Generated by Sollya with: |
| 100 | + // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 2^-2]); |
| 101 | + // The last coefficient 0x1.bd114ep-6f has been changed to 0x1.bd114ep-5f |
| 102 | + // for better accuracy. |
| 103 | + float p = fputil::polyeval(x_sq, 1.0f, -0x1.555552p-3f, 0x1.332f6ap-4f, |
| 104 | + -0x1.6c53dep-5f, 0x1.bd114ep-5f); |
| 105 | + |
| 106 | + return fputil::cast<float16>(xf * p); |
| 107 | + } |
| 108 | + |
| 109 | + // General case: asinh(x) = ln(x + sqrt(x^2 + 1)) |
| 110 | + float sqrt_term = fputil::sqrt<float>(fputil::multiply_add(xf, xf, 1.0f)); |
| 111 | + return fputil::cast<float16>( |
| 112 | + x_sign * log_eval(fputil::multiply_add(xf, x_sign, sqrt_term))); |
| 113 | +} |
| 114 | + |
| 115 | +} // namespace math |
| 116 | + |
| 117 | +} // namespace LIBC_NAMESPACE_DECL |
| 118 | + |
| 119 | +#endif // LIBC_TYPES_HAS_FLOAT16 |
| 120 | + |
| 121 | +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF16_H |
0 commit comments