17494Sjkh/*
27494Sjkh * Single-precision log(1+x) function.
37494Sjkh *
47494Sjkh * Copyright (c) 2022-2023, Arm Limited.
57494Sjkh * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
67494Sjkh */
77494Sjkh
87494Sjkh#include "poly_scalar_f32.h"
97494Sjkh#include "math_config.h"
107494Sjkh#include "pl_sig.h"
117494Sjkh#include "pl_test.h"
127494Sjkh
137494Sjkh#define Ln2 (0x1.62e43p-1f)
147494Sjkh#define SignMask (0x80000000)
157494Sjkh
167494Sjkh/* Biased exponent of the largest float m for which m^8 underflows.  */
177494Sjkh#define M8UFLOW_BOUND_BEXP 112
187494Sjkh/* Biased exponent of the largest float for which we just return x.  */
197494Sjkh#define TINY_BOUND_BEXP 103
207494Sjkh
217494Sjkh#define C(i) __log1pf_data.coeffs[i]
227494Sjkh
237494Sjkhstatic inline float
247494Sjkheval_poly (float m, uint32_t e)
257494Sjkh{
267494Sjkh#ifdef LOG1PF_2U5
277494Sjkh
287494Sjkh  /* 2.5 ulp variant. Approximate log(1+m) on [-0.25, 0.5] using
2965437Sphantom     slightly modified Estrin scheme (no x^0 term, and x term is just x).  */
30142665Sphantom  float p_12 = fmaf (m, C (1), C (0));
317494Sjkh  float p_34 = fmaf (m, C (3), C (2));
327494Sjkh  float p_56 = fmaf (m, C (5), C (4));
337494Sjkh  float p_78 = fmaf (m, C (7), C (6));
347494Sjkh
3579754Sdd  float m2 = m * m;
3659460Sphantom  float p_02 = fmaf (m2, p_12, m);
3759460Sphantom  float p_36 = fmaf (m2, p_56, p_34);
387494Sjkh  float p_79 = fmaf (m2, C (8), p_78);
3984306Sru
407494Sjkh  float m4 = m2 * m2;
4135548Sache  float p_06 = fmaf (m4, p_36, p_02);
427494Sjkh
4379754Sdd  if (unlikely (e < M8UFLOW_BOUND_BEXP))
4479754Sdd    return p_06;
457494Sjkh
467494Sjkh  float m8 = m4 * m4;
477494Sjkh  return fmaf (m8, p_79, p_06);
487494Sjkh
497494Sjkh#elif defined(LOG1PF_1U3)
507494Sjkh
5179754Sdd  /* 1.3 ulp variant. Approximate log(1+m) on [-0.25, 0.5] using Horner
527494Sjkh     scheme. Our polynomial approximation for log1p has the form
537494Sjkh     x + C1 * x^2 + C2 * x^3 + C3 * x^4 + ...
547494Sjkh     Hence approximation has the form m + m^2 * P(m)
5565437Sphantom       where P(x) = C1 + C2 * x + C3 * x^2 + ... .  */
5679754Sdd  return fmaf (m, m * horner_8_f32 (m, __log1pf_data.coeffs), m);
577494Sjkh
587494Sjkh#else
597494Sjkh#error No log1pf approximation exists with the requested precision. Options are 13 or 25.
607494Sjkh#endif
61142665Sphantom}
62142665Sphantom
63142665Sphantomstatic inline uint32_t
64142665Sphantombiased_exponent (uint32_t ix)
65142665Sphantom{
66147402Sru  return (ix & 0x7f800000) >> 23;
67142665Sphantom}
68142665Sphantom
69142665Sphantom/* log1pf approximation using polynomial on reduced interval. Worst-case error
70142665Sphantom   when using Estrin is roughly 2.02 ULP:
71142665Sphantom   log1pf(0x1.21e13ap-2) got 0x1.fe8028p-3 want 0x1.fe802cp-3.  */
72142665Sphantomfloat
73142665Sphantomlog1pf (float x)
747494Sjkh{
7551577Sphantom  uint32_t ix = asuint (x);
767494Sjkh  uint32_t ia = ix & ~SignMask;
777494Sjkh  uint32_t ia12 = ia >> 20;
787494Sjkh  uint32_t e = biased_exponent (ix);
797494Sjkh
807494Sjkh  /* Handle special cases first.  */
8135548Sache  if (unlikely (ia12 >= 0x7f8 || ix >= 0xbf800000 || ix == 0x80000000
82142665Sphantom		|| e <= TINY_BOUND_BEXP))
83    {
84      if (ix == 0xff800000)
85	{
86	  /* x == -Inf => log1pf(x) =  NaN.  */
87	  return NAN;
88	}
89      if ((ix == 0x7f800000 || e <= TINY_BOUND_BEXP) && ia12 <= 0x7f8)
90	{
91	  /* |x| < TinyBound => log1p(x)  =  x.
92	      x ==       Inf => log1pf(x) = Inf.  */
93	  return x;
94	}
95      if (ix == 0xbf800000)
96	{
97	  /* x == -1.0 => log1pf(x) = -Inf.  */
98	  return __math_divzerof (-1);
99	}
100      if (ia12 >= 0x7f8)
101	{
102	  /* x == +/-NaN => log1pf(x) = NaN.  */
103	  return __math_invalidf (asfloat (ia));
104	}
105      /* x <    -1.0 => log1pf(x) = NaN.  */
106      return __math_invalidf (x);
107    }
108
109  /* With x + 1 = t * 2^k (where t = m + 1 and k is chosen such that m
110			   is in [-0.25, 0.5]):
111     log1p(x) = log(t) + log(2^k) = log1p(m) + k*log(2).
112
113     We approximate log1p(m) with a polynomial, then scale by
114     k*log(2). Instead of doing this directly, we use an intermediate
115     scale factor s = 4*k*log(2) to ensure the scale is representable
116     as a normalised fp32 number.  */
117
118  if (ix <= 0x3f000000 || ia <= 0x3e800000)
119    {
120      /* If x is in [-0.25, 0.5] then we can shortcut all the logic
121	 below, as k = 0 and m = x.  All we need is to return the
122	 polynomial.  */
123      return eval_poly (x, e);
124    }
125
126  float m = x + 1.0f;
127
128  /* k is used scale the input. 0x3f400000 is chosen as we are trying to
129     reduce x to the range [-0.25, 0.5]. Inside this range, k is 0.
130     Outside this range, if k is reinterpreted as (NOT CONVERTED TO) float:
131	 let k = sign * 2^p      where sign = -1 if x < 0
132					       1 otherwise
133	 and p is a negative integer whose magnitude increases with the
134	 magnitude of x.  */
135  int k = (asuint (m) - 0x3f400000) & 0xff800000;
136
137  /* By using integer arithmetic, we obtain the necessary scaling by
138     subtracting the unbiased exponent of k from the exponent of x.  */
139  float m_scale = asfloat (asuint (x) - k);
140
141  /* Scale up to ensure that the scale factor is representable as normalised
142     fp32 number (s in [2**-126,2**26]), and scale m down accordingly.  */
143  float s = asfloat (asuint (4.0f) - k);
144  m_scale = m_scale + fmaf (0.25f, s, -1.0f);
145
146  float p = eval_poly (m_scale, biased_exponent (asuint (m_scale)));
147
148  /* The scale factor to be applied back at the end - by multiplying float(k)
149     by 2^-23 we get the unbiased exponent of k.  */
150  float scale_back = (float) k * 0x1.0p-23f;
151
152  /* Apply the scaling back.  */
153  return fmaf (scale_back, Ln2, p);
154}
155
156PL_SIG (S, F, 1, log1p, -0.9, 10.0)
157PL_TEST_ULP (log1pf, 1.52)
158PL_TEST_SYM_INTERVAL (log1pf, 0.0, 0x1p-23, 50000)
159PL_TEST_SYM_INTERVAL (log1pf, 0x1p-23, 0.001, 50000)
160PL_TEST_SYM_INTERVAL (log1pf, 0.001, 1.0, 50000)
161PL_TEST_SYM_INTERVAL (log1pf, 1.0, inf, 5000)
162