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