1/*
2 * Double-precision log(1+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 "poly_scalar_f64.h"
9#include "math_config.h"
10#include "pl_sig.h"
11#include "pl_test.h"
12
13#define Ln2Hi 0x1.62e42fefa3800p-1
14#define Ln2Lo 0x1.ef35793c76730p-45
15#define HfRt2Top 0x3fe6a09e /* top32(asuint64(sqrt(2)/2)).  */
16#define OneMHfRt2Top                                                           \
17  0x00095f62 /* top32(asuint64(1)) - top32(asuint64(sqrt(2)/2)).  */
18#define OneTop12 0x3ff
19#define BottomMask 0xffffffff
20#define OneMHfRt2 0x3fd2bec333018866
21#define Rt2MOne 0x3fda827999fcef32
22#define AbsMask 0x7fffffffffffffff
23#define ExpM63 0x3c00
24
25static inline double
26eval_poly (double f)
27{
28  double f2 = f * f;
29  double f4 = f2 * f2;
30  double f8 = f4 * f4;
31  return estrin_18_f64 (f, f2, f4, f8, f8 * f8, __log1p_data.coeffs);
32}
33
34/* log1p approximation using polynomial on reduced interval. Largest
35   observed errors are near the lower boundary of the region where k
36   is 0.
37   Maximum measured error: 1.75ULP.
38   log1p(-0x1.2e1aea97b3e5cp-2) got -0x1.65fb8659a2f9p-2
39			       want -0x1.65fb8659a2f92p-2.  */
40double
41log1p (double x)
42{
43  uint64_t ix = asuint64 (x);
44  uint64_t ia = ix & AbsMask;
45  uint32_t ia16 = ia >> 48;
46
47  /* Handle special cases first.  */
48  if (unlikely (ia16 >= 0x7ff0 || ix >= 0xbff0000000000000
49		|| ix == 0x8000000000000000))
50    {
51      if (ix == 0x8000000000000000 || ix == 0x7ff0000000000000)
52	{
53	  /* x ==  -0 => log1p(x) =  -0.
54	     x == Inf => log1p(x) = Inf.  */
55	  return x;
56	}
57      if (ix == 0xbff0000000000000)
58	{
59	  /* x == -1 => log1p(x) = -Inf.  */
60	  return __math_divzero (-1);
61	  ;
62	}
63      if (ia16 >= 0x7ff0)
64	{
65	  /* x == +/-NaN => log1p(x) = NaN.  */
66	  return __math_invalid (asdouble (ia));
67	}
68      /* x  <      -1 => log1p(x) =  NaN.
69	 x ==    -Inf => log1p(x) =  NaN.  */
70      return __math_invalid (x);
71    }
72
73  /* With x + 1 = t * 2^k (where t = f + 1 and k is chosen such that f
74			   is in [sqrt(2)/2, sqrt(2)]):
75     log1p(x) = k*log(2) + log1p(f).
76
77     f may not be representable exactly, so we need a correction term:
78     let m = round(1 + x), c = (1 + x) - m.
79     c << m: at very small x, log1p(x) ~ x, hence:
80     log(1+x) - log(m) ~ c/m.
81
82     We therefore calculate log1p(x) by k*log2 + log1p(f) + c/m.  */
83
84  uint64_t sign = ix & ~AbsMask;
85  if (ia <= OneMHfRt2 || (!sign && ia <= Rt2MOne))
86    {
87      if (unlikely (ia16 <= ExpM63))
88	{
89	  /* If exponent of x <= -63 then shortcut the polynomial and avoid
90	     underflow by just returning x, which is exactly rounded in this
91	     region.  */
92	  return x;
93	}
94      /* If x is in [sqrt(2)/2 - 1, sqrt(2) - 1] then we can shortcut all the
95	 logic below, as k = 0 and f = x and therefore representable exactly.
96	 All we need is to return the polynomial.  */
97      return fma (x, eval_poly (x) * x, x);
98    }
99
100  /* Obtain correctly scaled k by manipulation in the exponent.  */
101  double m = x + 1;
102  uint64_t mi = asuint64 (m);
103  uint32_t u = (mi >> 32) + OneMHfRt2Top;
104  int32_t k = (int32_t) (u >> 20) - OneTop12;
105
106  /* Correction term c/m.  */
107  double cm = (x - (m - 1)) / m;
108
109  /* Reduce x to f in [sqrt(2)/2, sqrt(2)].  */
110  uint32_t utop = (u & 0x000fffff) + HfRt2Top;
111  uint64_t u_red = ((uint64_t) utop << 32) | (mi & BottomMask);
112  double f = asdouble (u_red) - 1;
113
114  /* Approximate log1p(x) on the reduced input using a polynomial. Because
115     log1p(0)=0 we choose an approximation of the form:
116	x + C0*x^2 + C1*x^3 + C2x^4 + ...
117     Hence approximation has the form f + f^2 * P(f)
118	where P(x) = C0 + C1*x + C2x^2 + ...  */
119  double p = fma (f, eval_poly (f) * f, f);
120
121  double kd = k;
122  double y = fma (Ln2Lo, kd, cm);
123  return y + fma (Ln2Hi, kd, p);
124}
125
126PL_SIG (S, D, 1, log1p, -0.9, 10.0)
127PL_TEST_ULP (log1p, 1.26)
128PL_TEST_SYM_INTERVAL (log1p, 0.0, 0x1p-23, 50000)
129PL_TEST_SYM_INTERVAL (log1p, 0x1p-23, 0.001, 50000)
130PL_TEST_SYM_INTERVAL (log1p, 0.001, 1.0, 50000)
131PL_TEST_SYM_INTERVAL (log1p, 1.0, inf, 5000)
132