sysdeps/ieee754/ldbl-128ibm/e_expl.c

   1 /* Quad-precision floating point e^x.
   2    Copyright (C) 1999-2013 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_private.h>
  68 #include <sysdeps/ieee754/ldbl-128/t_expl.h>
  69
  70 static const long double C[] = {
  71 /* Smallest integer x for which e^x overflows.  */
  72 #define himark C[0]
  73  709.78271289338399678773454114191496482L,
  74
  75 /* Largest integer x for which e^x underflows.  */
  76 #define lomark C[1]
  77 -744.44007192138126231410729844608163411L,
  78
  79 /* 3x2^96 */
  80 #define THREEp96 C[2]
  81  59421121885698253195157962752.0L,
  82
  83 /* 3x2^103 */
  84 #define THREEp103 C[3]
  85  30423614405477505635920876929024.0L,
  86
  87 /* 3x2^111 */
  88 #define THREEp111 C[4]
  89  7788445287802241442795744493830144.0L,
  90
  91 /* 1/ln(2) */
  92 #define M_1_LN2 C[5]
  93  1.44269504088896340735992468100189204L,
  94
  95 /* first 93 bits of ln(2) */
  96 #define M_LN2_0 C[6]
  97  0.693147180559945309417232121457981864L,
  98
  99 /* ln2_0 - ln(2) */
 100 #define M_LN2_1 C[7]
 101 -1.94704509238074995158795957333327386E-31L,
 102
 103 /* very small number */
 104 #define TINY C[8]
 105  1.0e-308L,
 106
 107 /* 2^16383 */
 108 #define TWO1023 C[9]
 109  8.988465674311579538646525953945123668E+307L,
 110
 111 /* 256 */
 112 #define TWO8 C[10]
 113  256.0L,
 114
 115 /* 32768 */
 116 #define TWO15 C[11]
 117  32768.0L,
 118
 119 /* Chebyshev polynom coefficients for (exp(x)-1)/x */
 120 #define P1 C[12]
 121 #define P2 C[13]
 122 #define P3 C[14]
 123 #define P4 C[15]
 124 #define P5 C[16]
 125 #define P6 C[17]
 126  0.5L,
 127  1.66666666666666666666666666666666683E-01L,
 128  4.16666666666666666666654902320001674E-02L,
 129  8.33333333333333333333314659767198461E-03L,
 130  1.38888888889899438565058018857254025E-03L,
 131  1.98412698413981650382436541785404286E-04L,
 132 };
 133
 134 long double
 135 __ieee754_expl (long double x)
 136 {
 137   /* Check for usual case.  */
 138   if (isless (x, himark) && isgreater (x, lomark))
 139     {
 140       int tval1, tval2, unsafe, n_i, exponent2;
 141       long double x22, n, result, xl;
 142       union ibm_extended_long_double ex2_u, scale_u;
 143       fenv_t oldenv;
 144
 145       feholdexcept (&oldenv);
 146 #ifdef FE_TONEAREST
 147       fesetround (FE_TONEAREST);
 148 #endif
 149
 150       n = __roundl (x*M_1_LN2);
 151       x = x-n*M_LN2_0;
 152       xl = n*M_LN2_1;
 153
 154       tval1 = __roundl (x*TWO8);
 155       x -= __expl_table[T_EXPL_ARG1+2*tval1];
 156       xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
 157
 158       tval2 = __roundl (x*TWO15);
 159       x -= __expl_table[T_EXPL_ARG2+2*tval2];
 160       xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
 161
 162       x = x + xl;
 163
 164       /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]).  */
 165       ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
 166                 * __expl_table[T_EXPL_RES2 + tval2];
 167       n_i = (int)n;
 168       /* 'unsafe' is 1 iff n_1 != 0.  */
 169       unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
 170       ex2_u.ieee.exponent += n_i >> unsafe;
 171       /* Fortunately, there are no subnormal lowpart doubles in
 172          __expl_table, only normal values and zeros.
 173          But after scaling it can be subnormal.  */
 174       exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe);
 175       if (ex2_u.ieee.exponent2 == 0)
 176         /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */;
 177       else if (exponent2 > 0)
 178         ex2_u.ieee.exponent2 = exponent2;
 179       else if (exponent2 <= -54)
 180         {
 181           ex2_u.ieee.exponent2 = 0;
 182           ex2_u.ieee.mantissa2 = 0;
 183           ex2_u.ieee.mantissa3 = 0;
 184         }
 185       else
 186         {
 187           static const double
 188             two54 = 1.80143985094819840000e+16, /* 4350000000000000 */
 189             twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */
 190           ex2_u.dd[1] *= two54;
 191           ex2_u.ieee.exponent2 += n_i >> unsafe;
 192           ex2_u.dd[1] *= twom54;
 193         }
 194
 195       /* Compute scale = 2^n_1.  */
 196       scale_u.d = 1.0L;
 197       scale_u.ieee.exponent += n_i - (n_i >> unsafe);
 198
 199       /* Approximate e^x2 - 1, using a seventh-degree polynomial,
 200          with maximum error in [-2^-16-2^-53,2^-16+2^-53]
 201          less than 4.8e-39.  */
 202       x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
 203
 204       /* Return result.  */
 205       fesetenv (&oldenv);
 206
 207       result = x22 * ex2_u.d + ex2_u.d;
 208
 209       /* Now we can test whether the result is ultimate or if we are unsure.
 210          In the later case we should probably call a mpn based routine to give
 211          the ultimate result.
 212          Empirically, this routine is already ultimate in about 99.9986% of
 213          cases, the test below for the round to nearest case will be false
 214          in ~ 99.9963% of cases.
 215          Without proc2 routine maximum error which has been seen is
 216          0.5000262 ulp.
 217
 218           union ieee854_long_double ex3_u;
 219
 220           #ifdef FE_TONEAREST
 221             fesetround (FE_TONEAREST);
 222           #endif
 223           ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
 224           ex2_u.d = result;
 225           ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
 226                                  - ex2_u.ieee.exponent;
 227           n_i = abs (ex3_u.d);
 228           n_i = (n_i + 1) / 2;
 229           fesetenv (&oldenv);
 230           #ifdef FE_TONEAREST
 231           if (fegetround () == FE_TONEAREST)
 232             n_i -= 0x4000;
 233           #endif
 234           if (!n_i) {
 235             return __ieee754_expl_proc2 (origx);
 236           }
 237        */
 238       if (!unsafe)
 239         return result;
 240       else
 241         return result * scale_u.d;
 242     }
 243   /* Exceptional cases:  */
 244   else if (isless (x, himark))
 245     {
 246       if (__isinfl (x))
 247         /* e^-inf == 0, with no error.  */
 248         return 0;
 249       else
 250         /* Underflow */
 251         return TINY * TINY;
 252     }
 253   else
 254     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
 255     return TWO1023*x;
 256 }
 257 strong_alias (__ieee754_expl, __expl_finite)