sysdeps/ieee754/ldbl-128/e_expl.c

   1 /* Quad-precision floating point e^x.
   2    Copyright (C) 1999-2019 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
   5    Partly based on double-precision code
   6    by Geoffrey Keating <geoffk@ozemail.com.au>
   7
   8    The GNU C Library is free software; you can redistribute it and/or
   9    modify it under the terms of the GNU Lesser General Public
  10    License as published by the Free Software Foundation; either
  11    version 2.1 of the License, or (at your option) any later version.
  12
  13    The GNU C Library is distributed in the hope that it will be useful,
  14    but WITHOUT ANY WARRANTY; without even the implied warranty of
  15    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16    Lesser General Public License for more details.
  17
  18    You should have received a copy of the GNU Lesser General Public
  19    License along with the GNU C Library; if not, see
  20    <http://www.gnu.org/licenses/>.  */
  21
  22 /* The basic design here is from
  23    Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with
  24    Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991,
  25    pp. 410-423.
  26
  27    We work with number pairs where the first number is the high part and
  28    the second one is the low part. Arithmetic with the high part numbers must
  29    be exact, without any roundoff errors.
  30
  31    The input value, X, is written as
  32    X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x
  33        - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl
  34
  35    where:
  36    - n is an integer, 16384 >= n >= -16495;
  37    - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205
  38    - t1 is an integer, 89 >= t1 >= -89
  39    - t2 is an integer, 65 >= t2 >= -65
  40    - |arg1[t1]-t1/256.0| < 2^-53
  41    - |arg2[t2]-t2/32768.0| < 2^-53
  42    - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53
  43
  44    Then e^x is approximated as
  45
  46    e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
  47                + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
  48                  * p (x + xl + n * ln(2)_1))
  49    where:
  50    - p(x) is a polynomial approximating e(x)-1
  51    - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table
  52    - e^(arg2[t2]_0 + arg2[t2]_1) likewise
  53    - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1.
  54
  55    If it happens that n_1 == 0 (this is the usual case), that multiplication
  56    is omitted.
  57    */
  58
  59 #ifndef _GNU_SOURCE
  60 #define _GNU_SOURCE
  61 #endif
  62 #include <float.h>
  63 #include <ieee754.h>
  64 #include <math.h>
  65 #include <fenv.h>
  66 #include <inttypes.h>
  67 #include <math-barriers.h>
  68 #include <math_private.h>
  69 #include <math-underflow.h>
  70 #include <stdlib.h>
  71 #include "t_expl.h"
  72
  73 static const _Float128 C[] = {
  74 /* Smallest integer x for which e^x overflows.  */
  75 #define himark C[0]
  76  L(11356.523406294143949491931077970765),
  77
  78 /* Largest integer x for which e^x underflows.  */
  79 #define lomark C[1]
  80 L(-11433.4627433362978788372438434526231),
  81
  82 /* 3x2^96 */
  83 #define THREEp96 C[2]
  84  L(59421121885698253195157962752.0),
  85
  86 /* 3x2^103 */
  87 #define THREEp103 C[3]
  88  L(30423614405477505635920876929024.0),
  89
  90 /* 3x2^111 */
  91 #define THREEp111 C[4]
  92  L(7788445287802241442795744493830144.0),
  93
  94 /* 1/ln(2) */
  95 #define M_1_LN2 C[5]
  96  L(1.44269504088896340735992468100189204),
  97
  98 /* first 93 bits of ln(2) */
  99 #define M_LN2_0 C[6]
 100  L(0.693147180559945309417232121457981864),
 101
 102 /* ln2_0 - ln(2) */
 103 #define M_LN2_1 C[7]
 104 L(-1.94704509238074995158795957333327386E-31),
 105
 106 /* very small number */
 107 #define TINY C[8]
 108  L(1.0e-4900),
 109
 110 /* 2^16383 */
 111 #define TWO16383 C[9]
 112  L(5.94865747678615882542879663314003565E+4931),
 113
 114 /* 256 */
 115 #define TWO8 C[10]
 116  256,
 117
 118 /* 32768 */
 119 #define TWO15 C[11]
 120  32768,
 121
 122 /* Chebyshev polynom coefficients for (exp(x)-1)/x */
 123 #define P1 C[12]
 124 #define P2 C[13]
 125 #define P3 C[14]
 126 #define P4 C[15]
 127 #define P5 C[16]
 128 #define P6 C[17]
 129  L(0.5),
 130  L(1.66666666666666666666666666666666683E-01),
 131  L(4.16666666666666666666654902320001674E-02),
 132  L(8.33333333333333333333314659767198461E-03),
 133  L(1.38888888889899438565058018857254025E-03),
 134  L(1.98412698413981650382436541785404286E-04),
 135 };
 136
 137 _Float128
 138 __ieee754_expl (_Float128 x)
 139 {
 140   /* Check for usual case.  */
 141   if (isless (x, himark) && isgreater (x, lomark))
 142     {
 143       int tval1, tval2, unsafe, n_i;
 144       _Float128 x22, n, t, result, xl;
 145       union ieee854_long_double ex2_u, scale_u;
 146       fenv_t oldenv;
 147
 148       feholdexcept (&oldenv);
 149 #ifdef FE_TONEAREST
 150       fesetround (FE_TONEAREST);
 151 #endif
 152
 153       /* Calculate n.  */
 154       n = x * M_1_LN2 + THREEp111;
 155       n -= THREEp111;
 156       x = x - n * M_LN2_0;
 157       xl = n * M_LN2_1;
 158
 159       /* Calculate t/256.  */
 160       t = x + THREEp103;
 161       t -= THREEp103;
 162
 163       /* Compute tval1 = t.  */
 164       tval1 = (int) (t * TWO8);
 165
 166       x -= __expl_table[T_EXPL_ARG1+2*tval1];
 167       xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
 168
 169       /* Calculate t/32768.  */
 170       t = x + THREEp96;
 171       t -= THREEp96;
 172
 173       /* Compute tval2 = t.  */
 174       tval2 = (int) (t * TWO15);
 175
 176       x -= __expl_table[T_EXPL_ARG2+2*tval2];
 177       xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
 178
 179       x = x + xl;
 180
 181       /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]).  */
 182       ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
 183                 * __expl_table[T_EXPL_RES2 + tval2];
 184       n_i = (int)n;
 185       /* 'unsafe' is 1 iff n_1 != 0.  */
 186       unsafe = abs(n_i) >= 15000;
 187       ex2_u.ieee.exponent += n_i >> unsafe;
 188
 189       /* Compute scale = 2^n_1.  */
 190       scale_u.d = 1;
 191       scale_u.ieee.exponent += n_i - (n_i >> unsafe);
 192
 193       /* Approximate e^x2 - 1, using a seventh-degree polynomial,
 194          with maximum error in [-2^-16-2^-53,2^-16+2^-53]
 195          less than 4.8e-39.  */
 196       x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
 197       math_force_eval (x22);
 198
 199       /* Return result.  */
 200       fesetenv (&oldenv);
 201
 202       result = x22 * ex2_u.d + ex2_u.d;
 203
 204       /* Now we can test whether the result is ultimate or if we are unsure.
 205          In the later case we should probably call a mpn based routine to give
 206          the ultimate result.
 207          Empirically, this routine is already ultimate in about 99.9986% of
 208          cases, the test below for the round to nearest case will be false
 209          in ~ 99.9963% of cases.
 210          Without proc2 routine maximum error which has been seen is
 211          0.5000262 ulp.
 212
 213           union ieee854_long_double ex3_u;
 214
 215           #ifdef FE_TONEAREST
 216             fesetround (FE_TONEAREST);
 217           #endif
 218           ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
 219           ex2_u.d = result;
 220           ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
 221                                  - ex2_u.ieee.exponent;
 222           n_i = abs (ex3_u.d);
 223           n_i = (n_i + 1) / 2;
 224           fesetenv (&oldenv);
 225           #ifdef FE_TONEAREST
 226           if (fegetround () == FE_TONEAREST)
 227             n_i -= 0x4000;
 228           #endif
 229           if (!n_i) {
 230             return __ieee754_expl_proc2 (origx);
 231           }
 232        */
 233       if (!unsafe)
 234         return result;
 235       else
 236         {
 237           result *= scale_u.d;
 238           math_check_force_underflow_nonneg (result);
 239           return result;
 240         }
 241     }
 242   /* Exceptional cases:  */
 243   else if (isless (x, himark))
 244     {
 245       if (isinf (x))
 246         /* e^-inf == 0, with no error.  */
 247         return 0;
 248       else
 249         /* Underflow */
 250         return TINY * TINY;
 251     }
 252   else
 253     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
 254     return TWO16383*x;
 255 }
 256 strong_alias (__ieee754_expl, __expl_finite)