Back to index

glibc  2.9
e_expl.c
Go to the documentation of this file.
00001 /* Quad-precision floating point e^x.
00002    Copyright (C) 1999 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
00005    Partly based on double-precision code
00006    by Geoffrey Keating <geoffk@ozemail.com.au>
00007 
00008    The GNU C Library is free software; you can redistribute it and/or
00009    modify it under the terms of the GNU Lesser General Public
00010    License as published by the Free Software Foundation; either
00011    version 2.1 of the License, or (at your option) any later version.
00012 
00013    The GNU C Library is distributed in the hope that it will be useful,
00014    but WITHOUT ANY WARRANTY; without even the implied warranty of
00015    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00016    Lesser General Public License for more details.
00017 
00018    You should have received a copy of the GNU Lesser General Public
00019    License along with the GNU C Library; if not, write to the Free
00020    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00021    02111-1307 USA.  */
00022 
00023 /* The basic design here is from
00024    Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with
00025    Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991,
00026    pp. 410-423.
00027 
00028    We work with number pairs where the first number is the high part and
00029    the second one is the low part. Arithmetic with the high part numbers must
00030    be exact, without any roundoff errors.
00031 
00032    The input value, X, is written as
00033    X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x
00034        - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl
00035 
00036    where:
00037    - n is an integer, 16384 >= n >= -16495;
00038    - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205
00039    - t1 is an integer, 89 >= t1 >= -89
00040    - t2 is an integer, 65 >= t2 >= -65
00041    - |arg1[t1]-t1/256.0| < 2^-53
00042    - |arg2[t2]-t2/32768.0| < 2^-53
00043    - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53
00044 
00045    Then e^x is approximated as
00046 
00047    e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
00048               + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
00049                * p (x + xl + n * ln(2)_1))
00050    where:
00051    - p(x) is a polynomial approximating e(x)-1
00052    - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table
00053    - e^(arg2[t2]_0 + arg2[t2]_1) likewise
00054    - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1.
00055 
00056    If it happens that n_1 == 0 (this is the usual case), that multiplication
00057    is omitted.
00058    */
00059 
00060 #ifndef _GNU_SOURCE
00061 #define _GNU_SOURCE
00062 #endif
00063 #include <float.h>
00064 #include <ieee754.h>
00065 #include <math.h>
00066 #include <fenv.h>
00067 #include <inttypes.h>
00068 #include <math_private.h>
00069 #include <stdlib.h>
00070 #include "t_expl.h"
00071 
00072 static const long double C[] = {
00073 /* Smallest integer x for which e^x overflows.  */
00074 #define himark C[0]
00075  11356.523406294143949491931077970765L,
00076  
00077 /* Largest integer x for which e^x underflows.  */
00078 #define lomark C[1]
00079 -11433.4627433362978788372438434526231L,
00080 
00081 /* 3x2^96 */
00082 #define THREEp96 C[2]
00083  59421121885698253195157962752.0L,
00084 
00085 /* 3x2^103 */
00086 #define THREEp103 C[3]
00087  30423614405477505635920876929024.0L,
00088 
00089 /* 3x2^111 */
00090 #define THREEp111 C[4]
00091  7788445287802241442795744493830144.0L,
00092 
00093 /* 1/ln(2) */
00094 #define M_1_LN2 C[5]
00095  1.44269504088896340735992468100189204L,
00096 
00097 /* first 93 bits of ln(2) */
00098 #define M_LN2_0 C[6]
00099  0.693147180559945309417232121457981864L,
00100 
00101 /* ln2_0 - ln(2) */
00102 #define M_LN2_1 C[7]
00103 -1.94704509238074995158795957333327386E-31L,
00104 
00105 /* very small number */
00106 #define TINY C[8]
00107  1.0e-4900L,
00108 
00109 /* 2^16383 */
00110 #define TWO16383 C[9]
00111  5.94865747678615882542879663314003565E+4931L,
00112 
00113 /* 256 */
00114 #define TWO8 C[10]
00115  256.0L,
00116 
00117 /* 32768 */
00118 #define TWO15 C[11]
00119  32768.0L,
00120 
00121 /* Chebyshev polynom coeficients for (exp(x)-1)/x */
00122 #define P1 C[12]
00123 #define P2 C[13]
00124 #define P3 C[14]
00125 #define P4 C[15]
00126 #define P5 C[16]
00127 #define P6 C[17]
00128  0.5L,
00129  1.66666666666666666666666666666666683E-01L,
00130  4.16666666666666666666654902320001674E-02L,
00131  8.33333333333333333333314659767198461E-03L,
00132  1.38888888889899438565058018857254025E-03L,
00133  1.98412698413981650382436541785404286E-04L,
00134 };
00135 
00136 long double
00137 __ieee754_expl (long double x)
00138 {
00139   /* Check for usual case.  */
00140   if (isless (x, himark) && isgreater (x, lomark))
00141     {
00142       int tval1, tval2, unsafe, n_i;
00143       long double x22, n, t, result, xl;
00144       union ieee854_long_double ex2_u, scale_u;
00145       fenv_t oldenv;
00146 
00147       feholdexcept (&oldenv);
00148 #ifdef FE_TONEAREST
00149       fesetround (FE_TONEAREST);
00150 #endif
00151 
00152       /* Calculate n.  */
00153       n = x * M_1_LN2 + THREEp111;
00154       n -= THREEp111;
00155       x = x - n * M_LN2_0;
00156       xl = n * M_LN2_1;
00157 
00158       /* Calculate t/256.  */
00159       t = x + THREEp103;
00160       t -= THREEp103;
00161 
00162       /* Compute tval1 = t.  */
00163       tval1 = (int) (t * TWO8);
00164 
00165       x -= __expl_table[T_EXPL_ARG1+2*tval1];
00166       xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
00167 
00168       /* Calculate t/32768.  */
00169       t = x + THREEp96;
00170       t -= THREEp96;
00171 
00172       /* Compute tval2 = t.  */
00173       tval2 = (int) (t * TWO15);
00174 
00175       x -= __expl_table[T_EXPL_ARG2+2*tval2];
00176       xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
00177 
00178       x = x + xl;
00179 
00180       /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]).  */
00181       ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
00182               * __expl_table[T_EXPL_RES2 + tval2];
00183       n_i = (int)n;
00184       /* 'unsafe' is 1 iff n_1 != 0.  */
00185       unsafe = abs(n_i) >= -LDBL_MIN_EXP - 1;
00186       ex2_u.ieee.exponent += n_i >> unsafe;
00187 
00188       /* Compute scale = 2^n_1.  */
00189       scale_u.d = 1.0L;
00190       scale_u.ieee.exponent += n_i - (n_i >> unsafe);
00191 
00192       /* Approximate e^x2 - 1, using a seventh-degree polynomial,
00193         with maximum error in [-2^-16-2^-53,2^-16+2^-53]
00194         less than 4.8e-39.  */
00195       x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
00196 
00197       /* Return result.  */
00198       fesetenv (&oldenv);
00199 
00200       result = x22 * ex2_u.d + ex2_u.d;
00201 
00202       /* Now we can test whether the result is ultimate or if we are unsure.
00203         In the later case we should probably call a mpn based routine to give
00204         the ultimate result.
00205         Empirically, this routine is already ultimate in about 99.9986% of
00206         cases, the test below for the round to nearest case will be false
00207         in ~ 99.9963% of cases.
00208         Without proc2 routine maximum error which has been seen is
00209         0.5000262 ulp.
00210 
00211          union ieee854_long_double ex3_u;
00212 
00213          #ifdef FE_TONEAREST
00214            fesetround (FE_TONEAREST);
00215          #endif
00216          ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
00217          ex2_u.d = result;
00218          ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
00219                              - ex2_u.ieee.exponent;
00220          n_i = abs (ex3_u.d);
00221          n_i = (n_i + 1) / 2;
00222          fesetenv (&oldenv);
00223          #ifdef FE_TONEAREST
00224          if (fegetround () == FE_TONEAREST)
00225            n_i -= 0x4000;
00226          #endif
00227          if (!n_i) {
00228            return __ieee754_expl_proc2 (origx);
00229          }
00230        */
00231       if (!unsafe)
00232        return result;
00233       else
00234        return result * scale_u.d;
00235     }
00236   /* Exceptional cases:  */
00237   else if (isless (x, himark))
00238     {
00239       if (__isinfl (x))
00240        /* e^-inf == 0, with no error.  */
00241        return 0;
00242       else
00243        /* Underflow */
00244        return TINY * TINY;
00245     }
00246   else
00247     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
00248     return TWO16383*x;
00249 }