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/expm1_2u5.c
48375 views
1
/*

2
 * Double-precision e^x - 1 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 "poly_scalar_f64.h"

9
#include "math_config.h"

10
#include "test_sig.h"

11
#include "test_defs.h"

12


13
#define InvLn2 0x1.71547652b82fep0

14
#define Ln2hi 0x1.62e42fefa39efp-1

15
#define Ln2lo 0x1.abc9e3b39803fp-56

16
#define Shift 0x1.8p52

17
/* 0x1p-51, below which expm1(x) is within 2 ULP of x.  */

18
#define TinyBound 0x3cc0000000000000

19
/* Above which expm1(x) overflows.  */

20
#define BigBound 0x1.63108c75a1937p+9

21
/* Below which expm1(x) rounds to 1.  */

22
#define NegBound -0x1.740bf7c0d927dp+9

23
#define AbsMask 0x7fffffffffffffff

24


25
/* Approximation for exp(x) - 1 using polynomial on a reduced interval.

26
   The maximum error observed error is 2.17 ULP:

27
   expm1(0x1.63f90a866748dp-2) got 0x1.a9af56603878ap-2

28
			      want 0x1.a9af566038788p-2.  */

29
double

30
expm1 (double x)

31
{

32
  uint64_t ix = asuint64 (x);

33
  uint64_t ax = ix & AbsMask;

34


35
  /* Tiny, +Infinity.  */

36
  if (ax <= TinyBound || ix == 0x7ff0000000000000)

37
    return x;

38


39
  /* +/-NaN.  */

40
  if (ax > 0x7ff0000000000000)

41
    return __math_invalid (x);

42


43
  /* Result is too large to be represented as a double.  */

44
  if (x >= 0x1.63108c75a1937p+9)

45
    return __math_oflow (0);

46


47
  /* Result rounds to -1 in double precision.  */

48
  if (x <= NegBound)

49
    return -1;

50


51
  /* Reduce argument to smaller range:

52
     Let i = round(x / ln2)

53
     and f = x - i * ln2, then f is in [-ln2/2, ln2/2].

54
     exp(x) - 1 = 2^i * (expm1(f) + 1) - 1

55
     where 2^i is exact because i is an integer.  */

56
  double j = fma (InvLn2, x, Shift) - Shift;

57
  int64_t i = j;

58
  double f = fma (j, -Ln2hi, x);

59
  f = fma (j, -Ln2lo, f);

60


61
  /* Approximate expm1(f) using polynomial.

62
     Taylor expansion for expm1(x) has the form:

63
	 x + ax^2 + bx^3 + cx^4 ....

64
     So we calculate the polynomial P(f) = a + bf + cf^2 + ...

65
     and assemble the approximation expm1(f) ~= f + f^2 * P(f).  */

66
  double f2 = f * f;

67
  double f4 = f2 * f2;

68
  double p = fma (f2, estrin_10_f64 (f, f2, f4, f4 * f4, __expm1_poly), f);

69


70
  /* Assemble the result, using a slight rearrangement to achieve acceptable

71
     accuracy.

72
     expm1(x) ~= 2^i * (p + 1) - 1

73
     Let t = 2^(i - 1).  */

74
  double t = ldexp (0.5, i);

75
  /* expm1(x) ~= 2 * (p * t + (t - 1/2)).  */

76
  return 2 * fma (p, t, t - 0.5);

77
}

78


79
TEST_SIG (S, D, 1, expm1, -9.9, 9.9)

80
TEST_ULP (expm1, 1.68)

81
TEST_SYM_INTERVAL (expm1, 0, 0x1p-51, 1000)

82
TEST_INTERVAL (expm1, 0x1p-51, 0x1.63108c75a1937p+9, 100000)

83
TEST_INTERVAL (expm1, -0x1p-51, -0x1.740bf7c0d927dp+9, 100000)

84
TEST_INTERVAL (expm1, 0x1.63108c75a1937p+9, inf, 100)

85
TEST_INTERVAL (expm1, -0x1.740bf7c0d927dp+9, -inf, 100)

86


87