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 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 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 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