expf.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 28 /* 29 * ISO-IEC-10967-2: Elementary Numerical Functions 30 * Signature: 31 * double exp(double x) 32 * 33 * Spec: 34 * exp(1) = e 35 * exp(x) = 1 if x ∈ F and exp(x) ≠ eˣ 36 * exp(x) = 1 if x = -0 37 * exp(x) = +inf if x = +inf 38 * exp(x) = 0 if x = -inf 39 * exp(x) = eˣ 40 * 41 * exp(x) overflows if (approximately) x > ln(FLT_MAX) i.e., 88.72.. 42 */ 43 44 #include <stdint.h> 45 #include <libm_util_amd.h> 46 #include <libm_special.h> 47 48 #include <libm_macros.h> 49 #include <libm/types.h> 50 51 #include <libm/typehelper.h> 52 #include <libm/compiler.h> 53 54 55 #define EXPF_N 6 56 #define EXPF_POLY_DEGREE 4 57 #if EXPF_N == 5 58 #define EXPF_POLY_DEGREE 3 59 #elif EXPF_N == 4 60 #define EXPF_POLY_DEGREE 3 61 #endif 62 63 #define EXPF_TABLE_SIZE (1 << EXPF_N) 64 #define EXPF_MAX_POLY_DEGREE 4 65 66 #define EXP_Y_ZERO 2 67 #define EXP_Y_INF 3 68 /* 69 * expf_data.h needs following to be defined before include 70 * - EXPF_N 71 * - EXPF_POLY_DEGREE 72 * - EXPF_TABLE_SIZE 73 * - EXPF_MAX_POLY_DEGREE 74 */ 75 76 #include "expf_data.h" 77 78 static const struct expf_data expf_v2_data = { 79 .ln2by_tblsz = 0x1.62e42fefa39efp-7, 80 .tblsz_byln2 = 0x1.71547652b82fep+6, 81 .Huge = 0x1.8000000000000p+52, 82 #if EXPF_N == 6 83 //.table_v3 = __two_to_jby64, 84 #elif EXPF_N == 5 85 .table_v3 = &__two_to_jby32, 86 #endif 87 88 .poly = { 89 1.0, /* 1/1! = 1 */ 90 0x1.0000000000000p-1, /* 1/2! = 1/2 */ 91 0x1.5555555555555p-3, /* 1/3! = 1/6 */ 92 0x1.cacccaa4ba57cp-5, /* 1/4! = 1/24 */ 93 }, 94 }; 95 96 #define C1 expf_v2_data.poly[0] 97 #define C2 expf_v2_data.poly[1] 98 #define C3 expf_v2_data.poly[2] 99 #define C4 expf_v2_data.poly[3] 100 101 #define EXPF_LN2_BY_TBLSZ expf_v2_data.ln2by_tblsz 102 #define EXPF_TBLSZ_BY_LN2 expf_v2_data.tblsz_byln2 103 #define EXPF_HUGE expf_v2_data.Huge 104 #define EXPF_TABLE __two_to_jby64 105 106 #define EXPF_FARG_MIN -0x1.9fe368p6f /* log(0x1p-150) ~= -103.97 */ 107 #define EXPF_FARG_MAX 0x1.62e42ep6f /* log(0x1p128) ~= 88.72 */ 108 109 float _expf_special(float x, float y, uint32_t code); 110 111 static uint32_t 112 top12f(float x) 113 { 114 flt32_t f = {.f = x}; 115 return f.i >> 20; 116 } 117 118 /****************************************** 119 * Implementation Notes 120 * --------------------- 121 * 122 * 0. Choose 'N' as 5, EXPF_TBL_SZ = 2^N i.e 32 123 * 124 * 1. Argument Reduction 125 ******************************************/ 126 #undef EXPF_N 127 #define EXPF_N 6 128 129 #undef EXPF_TABLE_SIZE 130 #define EXPF_TABLE_SIZE (1 << EXPF_N) 131 132 float 133 ALM_PROTO_OPT(expf)(float x) 134 { 135 double_t q, dn, r, z; 136 uint64_t n, j; 137 138 uint32_t top = top12f(x); 139 140 if (unlikely (top > top12f(88.0f))) { 141 if(isnanf(x)) 142 return x; 143 144 if (asuint32(x) == asuint32(-INFINITY)) 145 return 0.0f; 146 147 if (x > EXPF_FARG_MAX){ 148 if(asuint32(x) == PINFBITPATT_SP32) 149 return asfloat(PINFBITPATT_SP32); 150 151 /* Raise FE_OVERFLOW, FE_INEXACT */ 152 return _expf_special(x, asfloat(PINFBITPATT_SP32), EXP_Y_INF); 153 } 154 155 if (x < EXPF_FARG_MIN){ 156 return _expf_special(x, 0.0, EXP_Y_ZERO);; 157 } 158 } 159 160 z = x * EXPF_TBLSZ_BY_LN2; 161 162 /* 163 * n = (int) scale(x) 164 * dn = (double) n 165 */ 166 #undef FAST_INTEGER_CONVERSION 167 #define FAST_INTEGER_CONVERSION 1 168 #if FAST_INTEGER_CONVERSION 169 dn = z + EXPF_HUGE; 170 171 n = asuint64(dn); 172 173 dn -= EXPF_HUGE; 174 #else 175 n = z; 176 dn = cast_i32_to_float(n); 177 178 #endif 179 180 r = x - dn * EXPF_LN2_BY_TBLSZ; 181 182 j = n % EXPF_TABLE_SIZE; 183 184 double_t qtmp = C2 + (C3 * r); 185 186 double_t r2 = r * r; 187 188 double_t tbl = asdouble(EXPF_TABLE[j] + (n << (52 - EXPF_N))); 189 190 q = r + (r2 * qtmp); 191 192 double_t result = tbl + tbl* q; 193 194 return (float_t)(result); 195 }