1 /* Exponential function. 2 Copyright (C) 2011-2021 Free Software Foundation, Inc. 3 4 This file is free software: you can redistribute it and/or modify 5 it under the terms of the GNU Lesser General Public License as 6 published by the Free Software Foundation; either version 3 of the 7 License, or (at your option) any later version. 8 9 This file is distributed in the hope that it will be useful, 10 but WITHOUT ANY WARRANTY; without even the implied warranty of 11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 GNU Lesser General Public License for more details. 13 14 You should have received a copy of the GNU Lesser General Public License 15 along with this program. If not, see <https://www.gnu.org/licenses/>. */ 16 17 #include <config.h> 18 19 /* Specification. */ 20 #include <math.h> 21 22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE 23 24 long double 25 expl (long double x) /* */ 26 { 27 return exp (x); 28 } 29 30 #else 31 32 # include <float.h> 33 34 /* gl_expl_table[i] = exp((i - 128) * log(2)/256). */ 35 extern const long double gl_expl_table[257]; 36 37 /* A value slightly larger than log(2). */ 38 #define LOG2_PLUS_EPSILON 0.6931471805599454L 39 40 /* Best possible approximation of log(2) as a 'long double'. */ 41 #define LOG2 0.693147180559945309417232121458176568075L 42 43 /* Best possible approximation of 1/log(2) as a 'long double'. */ 44 #define LOG2_INVERSE 1.44269504088896340735992468100189213743L 45 46 /* Best possible approximation of log(2)/256 as a 'long double'. */ 47 #define LOG2_BY_256 0.00270760617406228636491106297444600221904L 48 49 /* Best possible approximation of 256/log(2) as a 'long double'. */ 50 #define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L 51 52 /* The upper 32 bits of log(2)/256. */ 53 #define LOG2_BY_256_HI_PART 0.0027076061733168899081647396087646484375L 54 /* log(2)/256 - LOG2_HI_PART. */ 55 #define LOG2_BY_256_LO_PART \ 56 0.000000000000745396456746323365681353781544922399845L 57 58 long double 59 expl (long double x) /* */ 60 { 61 if (isnanl (x)) 62 return x; 63 64 if (x >= (long double) LDBL_MAX_EXP * LOG2_PLUS_EPSILON) 65 /* x > LDBL_MAX_EXP * log(2) 66 hence exp(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */ 67 return HUGE_VALL; 68 69 if (x <= (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * LOG2_PLUS_EPSILON) 70 /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) * log(2) 71 hence exp(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG), 72 underflows to zero. */ 73 return 0.0L; 74 75 /* Decompose x into 76 x = n * log(2) + m * log(2)/256 + y 77 where 78 n is an integer, 79 m is an integer, -128 <= m <= 128, 80 y is a number, |y| <= log(2)/512 + epsilon = 0.00135... 81 Then 82 exp(x) = 2^n * exp(m * log(2)/256) * exp(y) 83 The first factor is an ldexpl() call. 84 The second factor is a table lookup. 85 The third factor is computed 86 - either as sinh(y) + cosh(y) 87 where sinh(y) is computed through the power series: 88 sinh(y) = y + y^3/3! + y^5/5! + ... 89 and cosh(y) is computed as hypot(1, sinh(y)), 90 - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z)) 91 where z = y/2 92 and tanh(z) is computed through its power series: 93 tanh(z) = z 94 - 1/3 * z^3 95 + 2/15 * z^5 96 - 17/315 * z^7 97 + 62/2835 * z^9 98 - 1382/155925 * z^11 99 + 21844/6081075 * z^13 100 - 929569/638512875 * z^15 101 + ... 102 Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the 103 z^13 term is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we 104 can truncate the series after the z^11 term. 105 106 Given the usual bounds LDBL_MAX_EXP <= 16384, LDBL_MIN_EXP >= -16381, 107 LDBL_MANT_DIG <= 120, we can estimate x: -11440 <= x <= 11357. 108 This means, when dividing x by log(2), where we want x mod log(2) 109 to be precise to LDBL_MANT_DIG bits, we have to use an approximation 110 to log(2) that has 14+LDBL_MANT_DIG bits. */ 111 112 { 113 long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */ 114 /* n has at most 15 bits, nm therefore has at most 23 bits, therefore 115 n * LOG2_HI_PART is computed exactly, and n * LOG2_LO_PART is computed 116 with an absolute error < 2^15 * 2e-10 * 2^-LDBL_MANT_DIG. */ 117 long double y_tmp = x - nm * LOG2_BY_256_HI_PART; 118 long double y = y_tmp - nm * LOG2_BY_256_LO_PART; 119 long double z = 0.5L * y; 120 121 /* Coefficients of the power series for tanh(z). */ 122 #define TANH_COEFF_1 1.0L 123 #define TANH_COEFF_3 -0.333333333333333333333333333333333333334L 124 #define TANH_COEFF_5 0.133333333333333333333333333333333333334L 125 #define TANH_COEFF_7 -0.053968253968253968253968253968253968254L 126 #define TANH_COEFF_9 0.0218694885361552028218694885361552028218L 127 #define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L 128 #define TANH_COEFF_13 0.00359212803657248101692546136990581435026L 129 #define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L 130 131 long double z2 = z * z; 132 long double tanh_z = 133 (((((TANH_COEFF_11 134 * z2 + TANH_COEFF_9) 135 * z2 + TANH_COEFF_7) 136 * z2 + TANH_COEFF_5) 137 * z2 + TANH_COEFF_3) 138 * z2 + TANH_COEFF_1) 139 * z; 140 141 long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z); 142 143 int n = (int) roundl (nm * (1.0L / 256.0L)); 144 int m = (int) nm - 256 * n; 145 146 return ldexpl (gl_expl_table[128 + m] * exp_y, n); 147 } 148 } 149 150 #endif