1/*
2 * Single-precision sinh(x) function.
3 *
4 * Copyright (c) 2022-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8#include "math_config.h"
9#include "pl_sig.h"
10#include "pl_test.h"
11
12#define AbsMask 0x7fffffff
13#define Half 0x3f000000
14#define Expm1OFlowLimit                                                        \
15  0x42b17218 /* 0x1.62e43p+6, 2^7*ln2, minimum value for which expm1f          \
16		overflows.  */
17#define OFlowLimit                                                             \
18  0x42b2d4fd /* 0x1.65a9fap+6, minimum positive value for which sinhf should   \
19		overflow.  */
20
21float
22optr_aor_exp_f32 (float);
23
24/* Approximation for single-precision sinh(x) using expm1.
25   sinh(x) = (exp(x) - exp(-x)) / 2.
26   The maximum error is 2.26 ULP:
27   sinhf(0x1.e34a9ep-4) got 0x1.e469ep-4 want 0x1.e469e4p-4.  */
28float
29sinhf (float x)
30{
31  uint32_t ix = asuint (x);
32  uint32_t iax = ix & AbsMask;
33  float ax = asfloat (iax);
34  uint32_t sign = ix & ~AbsMask;
35  float halfsign = asfloat (Half | sign);
36
37  if (unlikely (iax >= Expm1OFlowLimit))
38    {
39      /* Special values and overflow.  */
40      if (iax >= 0x7fc00001 || iax == 0x7f800000)
41	return x;
42      if (iax >= 0x7f800000)
43	return __math_invalidf (x);
44      if (iax >= OFlowLimit)
45	return __math_oflowf (sign);
46
47      /* expm1f overflows a little before sinhf, (~88.7 vs ~89.4). We have to
48	 fill this gap by using a different algorithm, in this case we use a
49	 double-precision exp helper. For large x sinh(x) dominated by exp(x),
50	 however we cannot compute exp without overflow either. We use the
51	 identity:
52	 exp(a) = (exp(a / 2)) ^ 2.
53	 to compute sinh(x) ~= (exp(|x| / 2)) ^ 2 / 2    for x > 0
54			    ~= (exp(|x| / 2)) ^ 2 / -2   for x < 0.
55	 Greatest error in this region is 1.89 ULP:
56	 sinhf(0x1.65898cp+6) got 0x1.f00aep+127  want 0x1.f00adcp+127.  */
57      float e = optr_aor_exp_f32 (ax / 2);
58      return (e * halfsign) * e;
59    }
60
61  /* Use expm1f to retain acceptable precision for small numbers.
62     Let t = e^(|x|) - 1.  */
63  float t = expm1f (ax);
64  /* Then sinh(x) = (t + t / (t + 1)) / 2   for x > 0
65		    (t + t / (t + 1)) / -2  for x < 0.  */
66  return (t + t / (t + 1)) * halfsign;
67}
68
69PL_SIG (S, F, 1, sinh, -10.0, 10.0)
70PL_TEST_ULP (sinhf, 1.76)
71PL_TEST_SYM_INTERVAL (sinhf, 0, 0x1.62e43p+6, 100000)
72PL_TEST_SYM_INTERVAL (sinhf, 0x1.62e43p+6, 0x1.65a9fap+6, 100)
73PL_TEST_SYM_INTERVAL (sinhf, 0x1.65a9fap+6, inf, 100)
74