1 /*
2 * Single-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_f32.h"
9 #include "math_config.h"
10 #include "pl_sig.h"
11 #include "pl_test.h"
12
13 #define Ln2 (0x1.62e43p-1f)
14 #define SignMask (0x80000000)
15
16 /* Biased exponent of the largest float m for which m^8 underflows. */
17 #define M8UFLOW_BOUND_BEXP 112
18 /* Biased exponent of the largest float for which we just return x. */
19 #define TINY_BOUND_BEXP 103
20
21 #define C(i) __log1pf_data.coeffs[i]
22
23 static inline float
eval_poly(float m,uint32_t e)24 eval_poly (float m, uint32_t e)
25 {
26 #ifdef LOG1PF_2U5
27
28 /* 2.5 ulp variant. Approximate log(1+m) on [-0.25, 0.5] using
29 slightly modified Estrin scheme (no x^0 term, and x term is just x). */
30 float p_12 = fmaf (m, C (1), C (0));
31 float p_34 = fmaf (m, C (3), C (2));
32 float p_56 = fmaf (m, C (5), C (4));
33 float p_78 = fmaf (m, C (7), C (6));
34
35 float m2 = m * m;
36 float p_02 = fmaf (m2, p_12, m);
37 float p_36 = fmaf (m2, p_56, p_34);
38 float p_79 = fmaf (m2, C (8), p_78);
39
40 float m4 = m2 * m2;
41 float p_06 = fmaf (m4, p_36, p_02);
42
43 if (unlikely (e < M8UFLOW_BOUND_BEXP))
44 return p_06;
45
46 float m8 = m4 * m4;
47 return fmaf (m8, p_79, p_06);
48
49 #elif defined(LOG1PF_1U3)
50
51 /* 1.3 ulp variant. Approximate log(1+m) on [-0.25, 0.5] using Horner
52 scheme. Our polynomial approximation for log1p has the form
53 x + C1 * x^2 + C2 * x^3 + C3 * x^4 + ...
54 Hence approximation has the form m + m^2 * P(m)
55 where P(x) = C1 + C2 * x + C3 * x^2 + ... . */
56 return fmaf (m, m * horner_8_f32 (m, __log1pf_data.coeffs), m);
57
58 #else
59 #error No log1pf approximation exists with the requested precision. Options are 13 or 25.
60 #endif
61 }
62
63 static inline uint32_t
biased_exponent(uint32_t ix)64 biased_exponent (uint32_t ix)
65 {
66 return (ix & 0x7f800000) >> 23;
67 }
68
69 /* log1pf approximation using polynomial on reduced interval. Worst-case error
70 when using Estrin is roughly 2.02 ULP:
71 log1pf(0x1.21e13ap-2) got 0x1.fe8028p-3 want 0x1.fe802cp-3. */
72 float
log1pf(float x)73 log1pf (float x)
74 {
75 uint32_t ix = asuint (x);
76 uint32_t ia = ix & ~SignMask;
77 uint32_t ia12 = ia >> 20;
78 uint32_t e = biased_exponent (ix);
79
80 /* Handle special cases first. */
81 if (unlikely (ia12 >= 0x7f8 || ix >= 0xbf800000 || ix == 0x80000000
82 || 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
156 PL_SIG (S, F, 1, log1p, -0.9, 10.0)
157 PL_TEST_ULP (log1pf, 1.52)
158 PL_TEST_SYM_INTERVAL (log1pf, 0.0, 0x1p-23, 50000)
159 PL_TEST_SYM_INTERVAL (log1pf, 0x1p-23, 0.001, 50000)
160 PL_TEST_SYM_INTERVAL (log1pf, 0.001, 1.0, 50000)
161 PL_TEST_SYM_INTERVAL (log1pf, 1.0, inf, 5000)
162