Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
freebsd
GitHub Repository: freebsd/freebsd-src
 Path: blob/main/contrib/arm-optimized-routines/math/aarch64/experimental/atanh_3u.c
48378 views
1
/*

2
 * Double-precision atanh(x) function.

3
 *

4
 * Copyright (c) 2022-2024, Arm Limited.

5
 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception

6
 */

7


8
#include "math_config.h"

9
#include "poly_scalar_f64.h"

10
#include "test_sig.h"

11
#include "test_defs.h"

12


13
#define AbsMask 0x7fffffffffffffff

14
#define Half 0x3fe0000000000000

15
#define One 0x3ff0000000000000

16
#define Ln2Hi 0x1.62e42fefa3800p-1

17
#define Ln2Lo 0x1.ef35793c76730p-45

18
#define OneMHfRt2Top                                                          \

19
  0x00095f62 /* top32(asuint64(1)) - top32(asuint64(sqrt(2)/2)).  */

20
#define OneTop12 0x3ff

21
#define HfRt2Top 0x3fe6a09e /* top32(asuint64(sqrt(2)/2)).  */

22
#define BottomMask 0xffffffff

23


24
static inline double

25
log1p_inline (double x)

26
{

27
  /* Helper for calculating log(1 + x) using order-18 polynomial on a reduced

28
     interval. Copied from log1p_2u.c, with no special-case handling. See that

29
     file for details of the algorithm.  */

30
  double m = x + 1;

31
  uint64_t mi = asuint64 (m);

32


33
  /* Decompose x + 1 into (f + 1) * 2^k, with k chosen such that f is in

34
     [sqrt(2)/2, sqrt(2)].  */

35
  uint32_t u = (mi >> 32) + OneMHfRt2Top;

36
  int32_t k = (int32_t) (u >> 20) - OneTop12;

37
  uint32_t utop = (u & 0x000fffff) + HfRt2Top;

38
  uint64_t u_red = ((uint64_t) utop << 32) | (mi & BottomMask);

39
  double f = asdouble (u_red) - 1;

40


41
  /* Correction term for round-off in f.  */

42
  double cm = (x - (m - 1)) / m;

43


44
  /* Approximate log1p(f) with polynomial.  */

45
  double f2 = f * f;

46
  double f4 = f2 * f2;

47
  double f8 = f4 * f4;

48
  double p = fma (

49
      f, estrin_18_f64 (f, f2, f4, f8, f8 * f8, __log1p_data.coeffs) * f, f);

50


51
  /* Recombine log1p(x) = k*log2 + log1p(f) + c/m.  */

52
  double kd = k;

53
  double y = fma (Ln2Lo, kd, cm);

54
  return y + fma (Ln2Hi, kd, p);

55
}

56


57
/* Approximation for double-precision inverse tanh(x), using a simplified

58
   version of log1p. Greatest observed error is 3.00 ULP:

59
   atanh(0x1.e58f3c108d714p-4) got 0x1.e7da77672a647p-4

60
			      want 0x1.e7da77672a64ap-4.  */

61
double

62
atanh (double x)

63
{

64
  uint64_t ix = asuint64 (x);

65
  uint64_t sign = ix & ~AbsMask;

66
  uint64_t ia = ix & AbsMask;

67


68
  if (unlikely (ia == One))

69
    return __math_divzero (sign >> 32);

70


71
  if (unlikely (ia > One))

72
    return __math_invalid (x);

73


74
  double halfsign = asdouble (Half | sign);

75
  double ax = asdouble (ia);

76
  return halfsign * log1p_inline ((2 * ax) / (1 - ax));

77
}

78


79
TEST_SIG (S, D, 1, atanh, -1.0, 1.0)

80
TEST_ULP (atanh, 3.00)

81
TEST_SYM_INTERVAL (atanh, 0, 0x1p-23, 10000)

82
TEST_SYM_INTERVAL (atanh, 0x1p-23, 1, 90000)

83
TEST_SYM_INTERVAL (atanh, 1, inf, 100)

84


85