1//===-- Single-precision asinh 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-exception
6//
7//===----------------------------------------------------------------------===//
8
9#include "src/math/asinhf.h"
10#include "src/__support/FPUtil/FPBits.h"
11#include "src/__support/FPUtil/PolyEval.h"
12#include "src/__support/FPUtil/multiply_add.h"
13#include "src/__support/FPUtil/sqrt.h"
14#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
15#include "src/math/generic/common_constants.h"
16#include "src/math/generic/explogxf.h"
17
18namespace LIBC_NAMESPACE {
19
20LLVM_LIBC_FUNCTION(float, asinhf, (float x)) {
21 using FPBits_t = typename fputil::FPBits<float>;
22 FPBits_t xbits(x);
23 uint32_t x_u = xbits.uintval();
24 uint32_t x_abs = xbits.abs().uintval();
25
26 // |x| <= 2^-3
27 if (LIBC_UNLIKELY(x_abs <= 0x3e80'0000U)) {
28 // |x| <= 2^-26
29 if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
30 return static_cast<float>(LIBC_UNLIKELY(x_abs == 0)
31 ? x
32 : (x - 0x1.5555555555555p-3 * x * x * x));
33 }
34
35 double x_d = x;
36 double x_sq = x_d * x_d;
37 // Generated by Sollya with:
38 // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16], [|D...|],
39 // [0, 2^-2]);
40 double p = fputil::polyeval(
41 x: x_sq, a0: 0.0, a: -0x1.555555555551ep-3, a: 0x1.3333333325495p-4,
42 a: -0x1.6db6db5a7622bp-5, a: 0x1.f1c70f82928c6p-6, a: -0x1.6e893934266b7p-6,
43 a: 0x1.1c0b41d3fbe78p-6, a: -0x1.c0f47810b3c4fp-7, a: 0x1.2c8602690143dp-7);
44 return static_cast<float>(fputil::multiply_add(x: x_d, y: p, z: x_d));
45 }
46
47 const double SIGN[2] = {1.0, -1.0};
48 double x_sign = SIGN[x_u >> 31];
49 double x_d = x;
50
51 // Helper functions to set results for exceptional cases.
52 auto round_result_slightly_down = [x_sign](float r) -> float {
53 return fputil::multiply_add(x: static_cast<float>(x_sign), y: r,
54 z: static_cast<float>(x_sign) * (-0x1.0p-24f));
55 };
56 auto round_result_slightly_up = [x_sign](float r) -> float {
57 return fputil::multiply_add(x: static_cast<float>(x_sign), y: r,
58 z: static_cast<float>(x_sign) * 0x1.0p-24f);
59 };
60
61 if (LIBC_UNLIKELY(x_abs >= 0x4bdd'65a5U)) {
62 if (LIBC_UNLIKELY(xbits.is_inf_or_nan()))
63 return x;
64
65 // Exceptional cases when x > 2^24.
66 switch (x_abs) {
67 case 0x4bdd65a5: // |x| = 0x1.bacb4ap24f
68 return round_result_slightly_down(0x1.1e0696p4f);
69 case 0x4c803f2c: // |x| = 0x1.007e58p26f
70 return round_result_slightly_down(0x1.2b786cp4f);
71 case 0x4f8ffb03: // |x| = 0x1.1ff606p32f
72 return round_result_slightly_up(0x1.6fdd34p4f);
73 case 0x5c569e88: // |x| = 0x1.ad3d1p57f
74 return round_result_slightly_up(0x1.45c146p5f);
75 case 0x5e68984e: // |x| = 0x1.d1309cp61f
76 return round_result_slightly_up(0x1.5c9442p5f);
77 case 0x655890d3: // |x| = 0x1.b121a6p75f
78 return round_result_slightly_down(0x1.a9a3f2p5f);
79 case 0x65de7ca6: // |x| = 0x1.bcf94cp76f
80 return round_result_slightly_up(0x1.af66cp5f);
81 case 0x6eb1a8ec: // |x| = 0x1.6351d8p94f
82 return round_result_slightly_down(0x1.08b512p6f);
83 case 0x7997f30a: // |x| = 0x1.2fe614p116f
84 return round_result_slightly_up(0x1.451436p6f);
85 }
86 } else {
87 // Exceptional cases when x < 2^24.
88 if (LIBC_UNLIKELY(x_abs == 0x45abaf26)) {
89 // |x| = 0x1.575e4cp12f
90 return round_result_slightly_down(0x1.29becap3f);
91 }
92 if (LIBC_UNLIKELY(x_abs == 0x49d29048)) {
93 // |x| = 0x1.a5209p20f
94 return round_result_slightly_down(0x1.e1b92p3f);
95 }
96 }
97
98 // asinh(x) = log(x + sqrt(x^2 + 1))
99 return static_cast<float>(
100 x_sign *
101 log_eval(x: fputil::multiply_add(
102 x: x_d, y: x_sign, z: fputil::sqrt(x: fputil::multiply_add(x: x_d, y: x_d, z: 1.0)))));
103}
104
105} // namespace LIBC_NAMESPACE
106

source code of libc/src/math/generic/asinhf.cpp