lib/expl.c

   1 /* Exponential function.
   2    Copyright (C) 2011-2018 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program 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 General Public License for more details.
  13
  14    You should have received a copy of the GNU 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