| /* |
| * Copyright (c) 2014 Advanced Micro Devices, Inc. |
| * |
| * Permission is hereby granted, free of charge, to any person obtaining a copy |
| * of this software and associated documentation files (the "Software"), to deal |
| * in the Software without restriction, including without limitation the rights |
| * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
| * copies of the Software, and to permit persons to whom the Software is |
| * furnished to do so, subject to the following conditions: |
| * |
| * The above copyright notice and this permission notice shall be included in |
| * all copies or substantial portions of the Software. |
| * |
| * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
| * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
| * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
| * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
| * THE SOFTWARE. |
| */ |
| |
| #include "math64.h" |
| |
| // Algorithm: |
| // |
| // e^x = 2^(x/ln(2)) = 2^(x*(64/ln(2))/64) |
| // |
| // x*(64/ln(2)) = n + f, |f| <= 0.5, n is integer |
| // n = 64*m + j, 0 <= j < 64 |
| // |
| // e^x = 2^((64*m + j + f)/64) |
| // = (2^m) * (2^(j/64)) * 2^(f/64) |
| // = (2^m) * (2^(j/64)) * e^(f*(ln(2)/64)) |
| // |
| // f = x*(64/ln(2)) - n |
| // r = f*(ln(2)/64) = x - n*(ln(2)/64) |
| // |
| // e^x = (2^m) * (2^(j/64)) * e^r |
| // |
| // (2^(j/64)) is precomputed |
| // |
| // e^r = 1 + r + (r^2)/2! + (r^3)/3! + (r^4)/4! + (r^5)/5! |
| // e^r = 1 + q |
| // |
| // q = r + (r^2)/2! + (r^3)/3! + (r^4)/4! + (r^5)/5! |
| // |
| // e^x = (2^m) * ( (2^(j/64)) + q*(2^(j/64)) ) |
| |
| __attribute__((overloadable, always_inline, weak)) double |
| #if defined COMPILING_EXP2 |
| exp2(double x) |
| #elif defined COMPILING_EXP10 |
| exp10(double x) |
| #else |
| exp(double x) |
| #endif |
| { |
| USE_TABLE(double2, p_tbl, TWO_TO_JBY64_EP); |
| |
| #if defined(COMPILING_EXP2) |
| const double X_MAX = 1024.0; |
| const double X_MIN = -1074; |
| #elif defined(COMPILING_EXP10) |
| const double X_MAX = 0x1.34413509f79ffp+8; // 1024*ln(2)/ln(10) |
| const double X_MIN = -0x1.434e6420f4374p+8; // -1074*ln(2)/ln(10) |
| #else |
| const double X_MAX = 0x1.62e42fefa39efp+9; // 1024*ln(2) |
| const double X_MIN = -0x1.74910d52d3051p+9; // -1075*ln(2) |
| #endif |
| |
| #if defined(COMPILING_EXP2) |
| const double R_64 = 64.0; |
| const double R_1_BY_64 = 1.0 / 64.0; |
| const double R_LN2 = 0x1.62e42fefa39efp-1; // ln(2) |
| #elif defined(COMPILING_EXP10) |
| const double R_64_BY_LOG10_2 = 0x1.a934f0979a371p+7; // 64*ln(10)/ln(2) |
| const double R_LOG10_2_BY_64_LD = 0x1.3441350000000p-8; // head ln(2)/(64*ln(10)) |
| const double R_LOG10_2_BY_64_TL = 0x1.3ef3fde623e25p-37; // tail ln(2)/(64*ln(10)) |
| const double R_LN10 = 0x1.26bb1bbb55516p+1; // ln(10) |
| #else |
| const double R_64_BY_LOG2 = 0x1.71547652b82fep+6; // 64/ln(2) |
| const double R_LOG2_BY_64_LD = 0x1.62e42fefa0000p-7; // head ln(2)/64 |
| const double R_LOG2_BY_64_TL = 0x1.cf79abc9e3b39p-46; // tail ln(2)/64 |
| #endif |
| |
| #if defined(COMPILING_EXP2) |
| int n = convert_int(x * R_64); |
| #elif defined(COMPILING_EXP10) |
| int n = convert_int(x * R_64_BY_LOG10_2); |
| #else |
| int n = convert_int(x * R_64_BY_LOG2); |
| #endif |
| |
| double dn = (double)n; |
| |
| int j = n & 0x3f; |
| int m = n >> 6; |
| |
| #if defined(COMPILING_EXP2) |
| double r = R_LN2 * fma(-R_1_BY_64, dn, x); |
| #elif defined(COMPILING_EXP10) |
| double r = R_LN10 * fma(-R_LOG10_2_BY_64_TL, dn, fma(-R_LOG10_2_BY_64_LD, dn, x)); |
| #else |
| double r = fma(-R_LOG2_BY_64_TL, dn, fma(-R_LOG2_BY_64_LD, dn, x)); |
| #endif |
| |
| // 6 term tail of Taylor expansion of e^r |
| double z2 = r * fma(r, |
| fma(r, |
| fma(r, |
| fma(r, |
| fma(r, 0x1.6c16c16c16c17p-10, 0x1.1111111111111p-7), |
| 0x1.5555555555555p-5), |
| 0x1.5555555555555p-3), |
| 0x1.0000000000000p-1), |
| 1.0); |
| |
| double2 tv = p_tbl[j]; |
| z2 = fma(tv.s0 + tv.s1, z2, tv.s1) + tv.s0; |
| |
| int small_value = (m < -1022) || ((m == -1022) && (z2 < 1.0)); |
| |
| int n1 = m >> 2; |
| int n2 = m-n1; |
| double z3= z2 * as_double(((long)n1 + 1023) << 52); |
| z3 *= as_double(((long)n2 + 1023) << 52); |
| |
| z2 = ldexp(z2, m); |
| z2 = small_value ? z3: z2; |
| |
| z2 = isnan(x) ? x : z2; |
| |
| z2 = x > X_MAX ? as_double(PINFBITPATT_DP64) : z2; |
| z2 = x < X_MIN ? 0.0 : z2; |
| |
| return z2; |
| } |
| |