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

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