1/* 2 * Copyright (c) 2014 Advanced Micro Devices, Inc. 3 * 4 * Permission is hereby granted, free of charge, to any person obtaining a copy 5 * of this software and associated documentation files (the "Software"), to deal 6 * in the Software without restriction, including without limitation the rights 7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 8 * copies of the Software, and to permit persons to whom the Software is 9 * furnished to do so, subject to the following conditions: 10 * 11 * The above copyright notice and this permission notice shall be included in 12 * all copies or substantial portions of the Software. 13 * 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 20 * THE SOFTWARE. 21 */ 22 23#include <clc/clc.h> 24#include <clc/clcmacro.h> 25#include <clc/math/clc_mad.h> 26#include <clc/math/clc_subnormal_config.h> 27#include <clc/math/math.h> 28#include <clc/math/tables.h> 29#include <clc/relational/clc_isnan.h> 30 31// Algorithm: 32// 33// e^x = 2^(x/ln(2)) = 2^(x*(64/ln(2))/64) 34// 35// x*(64/ln(2)) = n + f, |f| <= 0.5, n is integer 36// n = 64*m + j, 0 <= j < 64 37// 38// e^x = 2^((64*m + j + f)/64) 39// = (2^m) * (2^(j/64)) * 2^(f/64) 40// = (2^m) * (2^(j/64)) * e^(f*(ln(2)/64)) 41// 42// f = x*(64/ln(2)) - n 43// r = f*(ln(2)/64) = x - n*(ln(2)/64) 44// 45// e^x = (2^m) * (2^(j/64)) * e^r 46// 47// (2^(j/64)) is precomputed 48// 49// e^r = 1 + r + (r^2)/2! + (r^3)/3! + (r^4)/4! + (r^5)/5! 50// e^r = 1 + q 51// 52// q = r + (r^2)/2! + (r^3)/3! + (r^4)/4! + (r^5)/5! 53// 54// e^x = (2^m) * ( (2^(j/64)) + q*(2^(j/64)) ) 55 56_CLC_DEF _CLC_OVERLOAD float __clc_exp10(float x) { 57 // 128*log2/log10 : 38.53183944498959 58 const float X_MAX = 0x1.344134p+5f; 59 // -149*log2/log10 : -44.8534693539332 60 const float X_MIN = -0x1.66d3e8p+5f; 61 // 64*log10/log2 : 212.6033980727912 62 const float R_64_BY_LOG10_2 = 0x1.a934f0p+7f; 63 // log2/(64 * log10) lead : 0.004699707 64 const float R_LOG10_2_BY_64_LD = 0x1.340000p-8f; 65 // log2/(64 * log10) tail : 0.00000388665057 66 const float R_LOG10_2_BY_64_TL = 0x1.04d426p-18f; 67 const float R_LN10 = 0x1.26bb1cp+1f; 68 69 int return_nan = __clc_isnan(x); 70 int return_inf = x > X_MAX; 71 int return_zero = x < X_MIN; 72 73 int n = convert_int(x * R_64_BY_LOG10_2); 74 75 float fn = (float)n; 76 int j = n & 0x3f; 77 int m = n >> 6; 78 int m2 = m << EXPSHIFTBITS_SP32; 79 float r; 80 81 r = R_LN10 * 82 __clc_mad(fn, -R_LOG10_2_BY_64_TL, __clc_mad(fn, -R_LOG10_2_BY_64_LD, x)); 83 84 // Truncated Taylor series for e^r 85 float z2 = __clc_mad(__clc_mad(__clc_mad(r, 0x1.555556p-5f, 0x1.555556p-3f), 86 r, 0x1.000000p-1f), 87 r * r, r); 88 89 float two_to_jby64 = USE_TABLE(exp_tbl, j); 90 z2 = __clc_mad(two_to_jby64, z2, two_to_jby64); 91 92 float z2s = z2 * as_float(0x1 << (m + 149)); 93 float z2n = as_float(as_int(z2) + m2); 94 z2 = m <= -126 ? z2s : z2n; 95 96 z2 = return_inf ? as_float(PINFBITPATT_SP32) : z2; 97 z2 = return_zero ? 0.0f : z2; 98 z2 = return_nan ? x : z2; 99 return z2; 100} 101_CLC_UNARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_exp10, float) 102 103#ifdef cl_khr_fp64 104_CLC_DEF _CLC_OVERLOAD double __clc_exp10(double x) { 105 // 1024*ln(2)/ln(10) 106 const double X_MAX = 0x1.34413509f79ffp+8; 107 // -1074*ln(2)/ln(10) 108 const double X_MIN = -0x1.434e6420f4374p+8; 109 // 64*ln(10)/ln(2) 110 const double R_64_BY_LOG10_2 = 0x1.a934f0979a371p+7; 111 // head ln(2)/(64*ln(10)) 112 const double R_LOG10_2_BY_64_LD = 0x1.3441350000000p-8; 113 // tail ln(2)/(64*ln(10)) 114 const double R_LOG10_2_BY_64_TL = 0x1.3ef3fde623e25p-37; 115 // ln(10) 116 const double R_LN10 = 0x1.26bb1bbb55516p+1; 117 118 int n = convert_int(x * R_64_BY_LOG10_2); 119 120 double dn = (double)n; 121 122 int j = n & 0x3f; 123 int m = n >> 6; 124 125 double r = 126 R_LN10 * fma(-R_LOG10_2_BY_64_TL, dn, fma(-R_LOG10_2_BY_64_LD, dn, x)); 127 128 // 6 term tail of Taylor expansion of e^r 129 double z2 = 130 r * 131 fma(r, 132 fma(r, 133 fma(r, 134 fma(r, fma(r, 0x1.6c16c16c16c17p-10, 0x1.1111111111111p-7), 135 0x1.5555555555555p-5), 136 0x1.5555555555555p-3), 137 0x1.0000000000000p-1), 138 1.0); 139 140 double2 tv = USE_TABLE(two_to_jby64_ep_tbl, j); 141 z2 = fma(tv.s0 + tv.s1, z2, tv.s1) + tv.s0; 142 143 int small_value = (m < -1022) || ((m == -1022) && (z2 < 1.0)); 144 145 int n1 = m >> 2; 146 int n2 = m - n1; 147 double z3 = z2 * as_double(((long)n1 + 1023) << 52); 148 z3 *= as_double(((long)n2 + 1023) << 52); 149 150 z2 = ldexp(z2, m); 151 z2 = small_value ? z3 : z2; 152 153 z2 = __clc_isnan(x) ? x : z2; 154 155 z2 = x > X_MAX ? as_double(PINFBITPATT_DP64) : z2; 156 z2 = x < X_MIN ? 0.0 : z2; 157 158 return z2; 159} 160_CLC_UNARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, double, __clc_exp10, double) 161#endif 162