Back to index

glibc  2.9
s_expm1l.c
Go to the documentation of this file.
00001 /*                                               expm1l.c
00002  *
00003  *     Exponential function, minus 1
00004  *      128-bit long double precision
00005  *
00006  *
00007  *
00008  * SYNOPSIS:
00009  *
00010  * long double x, y, expm1l();
00011  *
00012  * y = expm1l( x );
00013  *
00014  *
00015  *
00016  * DESCRIPTION:
00017  *
00018  * Returns e (2.71828...) raised to the x power, minus one.
00019  *
00020  * Range reduction is accomplished by separating the argument
00021  * into an integer k and fraction f such that
00022  *
00023  *     x    k  f
00024  *    e  = 2  e.
00025  *
00026  * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
00027  * in the basic range [-0.5 ln 2, 0.5 ln 2].
00028  *
00029  *
00030  * ACCURACY:
00031  *
00032  *                      Relative error:
00033  * arithmetic   domain     # trials      peak         rms
00034  *    IEEE    -79,+MAXLOG    100,000     1.7e-34     4.5e-35
00035  *
00036  */
00037 
00038 /* Copyright 2001 by Stephen L. Moshier 
00039 
00040     This library is free software; you can redistribute it and/or
00041     modify it under the terms of the GNU Lesser General Public
00042     License as published by the Free Software Foundation; either
00043     version 2.1 of the License, or (at your option) any later version.
00044 
00045     This library is distributed in the hope that it will be useful,
00046     but WITHOUT ANY WARRANTY; without even the implied warranty of
00047     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00048     Lesser General Public License for more details.
00049 
00050     You should have received a copy of the GNU Lesser General Public
00051     License along with this library; if not, write to the Free Software
00052     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00053 
00054 
00055 
00056 #include "math.h"
00057 #include "math_private.h"
00058 
00059 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
00060    -.5 ln 2  <  x  <  .5 ln 2
00061    Theoretical peak relative error = 8.1e-36  */
00062 
00063 static const long double
00064   P0 = 2.943520915569954073888921213330863757240E8L,
00065   P1 = -5.722847283900608941516165725053359168840E7L,
00066   P2 = 8.944630806357575461578107295909719817253E6L,
00067   P3 = -7.212432713558031519943281748462837065308E5L,
00068   P4 = 4.578962475841642634225390068461943438441E4L,
00069   P5 = -1.716772506388927649032068540558788106762E3L,
00070   P6 = 4.401308817383362136048032038528753151144E1L,
00071   P7 = -4.888737542888633647784737721812546636240E-1L,
00072   Q0 = 1.766112549341972444333352727998584753865E9L,
00073   Q1 = -7.848989743695296475743081255027098295771E8L,
00074   Q2 = 1.615869009634292424463780387327037251069E8L,
00075   Q3 = -2.019684072836541751428967854947019415698E7L,
00076   Q4 = 1.682912729190313538934190635536631941751E6L,
00077   Q5 = -9.615511549171441430850103489315371768998E4L,
00078   Q6 = 3.697714952261803935521187272204485251835E3L,
00079   Q7 = -8.802340681794263968892934703309274564037E1L,
00080   /* Q8 = 1.000000000000000000000000000000000000000E0 */
00081 /* C1 + C2 = ln 2 */
00082 
00083   C1 = 6.93145751953125E-1L,
00084   C2 = 1.428606820309417232121458176568075500134E-6L,
00085 /* ln (2^16384 * (1 - 2^-113)) */
00086   maxlog = 1.1356523406294143949491931077970764891253E4L,
00087 /* ln 2^-114 */
00088   minarg = -7.9018778583833765273564461846232128760607E1L, big = 2e4932L;
00089 
00090 
00091 long double
00092 __expm1l (long double x)
00093 {
00094   long double px, qx, xx;
00095   int32_t ix, sign;
00096   ieee854_long_double_shape_type u;
00097   int k;
00098 
00099   /* Detect infinity and NaN.  */
00100   u.value = x;
00101   ix = u.parts32.w0;
00102   sign = ix & 0x80000000;
00103   ix &= 0x7fffffff;
00104   if (ix >= 0x7fff0000)
00105     {
00106       /* Infinity. */
00107       if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
00108        {
00109          if (sign)
00110            return -1.0L;
00111          else
00112            return x;
00113        }
00114       /* NaN. No invalid exception. */
00115       return x;
00116     }
00117 
00118   /* expm1(+- 0) = +- 0.  */
00119   if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
00120     return x;
00121 
00122   /* Overflow.  */
00123   if (x > maxlog)
00124     return (big * big);
00125 
00126   /* Minimum value.  */
00127   if (x < minarg)
00128     return (4.0/big - 1.0L);
00129 
00130   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
00131   xx = C1 + C2;                    /* ln 2. */
00132   px = __floorl (0.5 + x / xx);
00133   k = px;
00134   /* remainder times ln 2 */
00135   x -= px * C1;
00136   x -= px * C2;
00137 
00138   /* Approximate exp(remainder ln 2).  */
00139   px = (((((((P7 * x
00140              + P6) * x
00141             + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
00142 
00143   qx = (((((((x
00144              + Q7) * x
00145             + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
00146 
00147   xx = x * x;
00148   qx = x + (0.5 * xx + xx * px / qx);
00149 
00150   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
00151 
00152   We have qx = exp(remainder ln 2) - 1, so
00153   exp(x) - 1 = 2^k (qx + 1) - 1
00154              = 2^k qx + 2^k - 1.  */
00155 
00156   px = __ldexpl (1.0L, k);
00157   x = px * qx + (px - 1.0);
00158   return x;
00159 }
00160 libm_hidden_def (__expm1l)
00161 weak_alias (__expm1l, expm1l)