logf.c
1 /* 2 * Copyright (C) 2008-2020 Advanced Micro Devices, Inc. All rights reserved. 3 * 4 * Redistribution and use in source and binary forms, with or without modification, 5 * are permitted provided that the following conditions are met: 6 * 1. Redistributions of source code must retain the above copyright notice, 7 * this list of conditions and the following disclaimer. 8 * 2. Redistributions in binary form must reproduce the above copyright notice, 9 * this list of conditions and the following disclaimer in the documentation 10 * and/or other materials provided with the distribution. 11 * 3. Neither the name of the copyright holder nor the names of its contributors 12 * may be used to endorse or promote products derived from this software without 13 * specific prior written permission. 14 * 15 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 17 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 18 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, 19 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 20 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 21 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 22 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 23 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 24 * POSSIBILITY OF SUCH DAMAGE. 25 * 26 */ 27 #include <stdint.h> 28 29 #include "libm_macros.h" 30 #include "libm_util_amd.h" 31 #include "libm_special.h" 32 33 #include <libm/types.h> 34 #include <libm/typehelper.h> 35 #include <libm/compiler.h> 36 37 #define MASK_MANT_ALL7 0x007f8000 38 #define MASK_MANT8 0x00008000 39 #define MASK_SIGN 0x7FFFFFFF 40 #define LOGF_N 8 41 #define LOGF_POLY_DEGREE 2 42 43 #define FLAG_X_ZERO 0x1 44 #define FLAG_X_NEG 0x2 45 #define FLAG_X_NAN 0x3 46 47 #include "logf_data.h" 48 49 extern float_t _logf_special(float_t x, float_t y, uint32_t code); 50 51 extern struct logf_table logf_lookup[]; 52 53 static struct { 54 float_t log2_head, log2_tail; 55 float_t poly[LOGF_POLY_DEGREE]; 56 //const float_t *tab; 57 float_t poly1[LOGF_POLY_DEGREE]; 58 } logf_data = { 59 .log2_head = 0x1.62e3p-1, 60 .log2_tail = 0x1.2fefa2p-17, 61 62 /* Polynomial constants */ 63 .poly = { 64 0.5f, /* C1= 1/2 */ 65 3.333333432674407958984375E-1f, /* C2 = 1/3 */ 66 }, 67 //.tab = logf_lookup, 68 69 /* Polynomial constants for cases near to 1 */ 70 .poly1 = { 71 8.33333333333317923934e-02f, /* A1 */ 72 1.25000000037717509602e-02, /* A2 */ 73 }, 74 }; 75 #define LOG2_HEAD logf_data.log2_head 76 #define LOG2_TAIL logf_data.log2_tail 77 78 static inline float_t 79 inline_log1pf(float_t x) 80 { 81 float_t r, r2, w; 82 83 r = x - 1.0f; 84 85 w = r / (2.0f + r); 86 87 float_t correction = w * r; 88 89 w += w; 90 91 float_t w2 = w * w; 92 93 #define A1 logf_data.poly1[0] 94 #define A2 logf_data.poly1[1] 95 96 r2 = (w * w2 * (A1 + A2 * w2)) - correction; 97 98 float_t f = r2 + r; 99 100 return f; 101 } 102 103 104 /* 105 * ISO-IEC-10967-2: Elementary Numerical Functions 106 * Signature: 107 * float logf(float x) 108 * 109 * Spec: 110 * logf(x) 111 * = logf(x) if x ∈ F and x > 0 112 * = x if x = qNaN 113 * = 0 if x = 1 114 * = -inf if x = (-0, 0} 115 * = NaN otherwise 116 * 117 * Implementation Notes: 118 * x = 2^m * (1 + A) where 1 <= A < 2 119 * 120 * x = 2^m * (1 + (G + g)) where 1 <= 1+G < 2 and g < 2^(-8) 121 * 122 * Let F = (1 + G) and f = g 123 * x = 2^m*(F + f), so 1 <= F < 2, f < 2^(-8) 124 * 125 * x = (2^m) * 2 * (F/2+f/2) 126 * So, 0.5 <= F/2 < 1, f/2 < 2^-9 127 * 128 * logf(x) = logf(2^m) + logf(2) + logf(F/2 + f/2) 129 * => A + B + C 130 * 131 * A = logf(2^m) = m*logf(2) .. (1) 132 * B = logf(2) .. (2) => Constant can be precalculated 133 * C = logf(F/2*(1 + f/F)) .. (3) 134 * => logf(F/2*(1 + f/F)) 135 * => logf(F/2) + logf(1 + f/F) 136 * 137 * logf(x) = logf(2^m) + logf(2) + logf(F/2) + logf(1 + f/F) 138 * 139 * (from log(a) + log(b) <=> log(a*b)) 140 * 141 * => m*logf(2) + logf(F) + logf(1 + r) where r = 1+f/F 142 */ 143 144 float 145 ALM_PROTO_OPT(logf)(float x) 146 { 147 uint32_t ux = asuint32(x); 148 149 if (unlikely (ux - 0x00800000 >= 0x7f800000 - 0x00800000)) 150 { 151 /* x < 0x1p-126 or inf or nan. */ 152 if (ux * 2 == 0) /* log(0) = -inf */ 153 return -1.0/0.0; 154 if (ux == 0x7f800000) /* log(inf) = inf */ 155 return x; 156 if ((ux & 0x80000000) || ux * 2 >= 0xff000000) 157 return sqrt(x); /* Return NaN */ 158 159 /* 160 * 'x' has to be denormal, Normalize it 161 * there is a possibility that only the last (23'rd) bit is set 162 * Hence multiply by 2^23 to bring to 1.xxxxx format. 163 */ 164 ux = asuint32(x * 0x1p23f); 165 ux -= 23 << 23; 166 } 167 168 int32_t expo = (ux >> EXPSHIFTBITS_SP32) - EMAX_SP32; 169 float_t f_expo = (float_t)expo; 170 171 #define NEAR_ONE_LO asuint32(1 - 0x1.0p-4) 172 #define NEAR_ONE_HI asuint32(1 + 0x1.1p-4) 173 174 /* Values very close to 1, e^(-1/16) <= x <= e^(1/16)*/ 175 if (unlikely(ux - NEAR_ONE_LO < NEAR_ONE_HI - NEAR_ONE_LO)) { 176 return inline_log1pf(x); 177 } 178 179 /* 180 * Here onwards, 'x' is neither -ve, nor close to 1 181 */ 182 uint32_t mant, mant1, idx; 183 float_t y, f, finv, r, r2, q, w; 184 185 mant = ux & MANTBITS_SP32; 186 mant1 = ux & MASK_MANT_ALL7; 187 /*This step is needed for more accuracy */ 188 /* mant1 += ((ux & MASK_MANT8) << 1); */ 189 190 idx = mant1 >> (EXPSHIFTBITS_SP32 - LOGF_N); 191 192 y = asfloat(mant |= HALFEXPBITS_SP32); 193 f = asfloat(mant1 |= HALFEXPBITS_SP32); 194 195 struct logf_table *tbl = &logf_lookup[idx]; 196 197 finv = asfloat(tbl->f_inv); 198 199 r = (f - y) * finv; 200 201 #define C1 logf_data.poly[0] 202 #define C2 logf_data.poly[1] 203 204 r2 = r * r; /* r^2 */ 205 q = r + (r2 * (C1 + r * C2)); 206 207 q = (f_expo * LOG2_TAIL) - q; 208 209 w = (LOG2_HEAD * f_expo) + asfloat(tbl->f_128_head); 210 211 q = (asfloat(tbl->f_128_tail) + q) + w; 212 213 return q; 214 } 215