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 #include "math.h"
00055 #include "math_private.h"
00056 #include <math_ldbl_opt.h>
00057 
00058 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
00059    -.5 ln 2  <  x  <  .5 ln 2
00060    Theoretical peak relative error = 8.1e-36  */
00061 
00062 static const long double
00063   P0 = 2.943520915569954073888921213330863757240E8L,
00064   P1 = -5.722847283900608941516165725053359168840E7L,
00065   P2 = 8.944630806357575461578107295909719817253E6L,
00066   P3 = -7.212432713558031519943281748462837065308E5L,
00067   P4 = 4.578962475841642634225390068461943438441E4L,
00068   P5 = -1.716772506388927649032068540558788106762E3L,
00069   P6 = 4.401308817383362136048032038528753151144E1L,
00070   P7 = -4.888737542888633647784737721812546636240E-1L,
00071   Q0 = 1.766112549341972444333352727998584753865E9L,
00072   Q1 = -7.848989743695296475743081255027098295771E8L,
00073   Q2 = 1.615869009634292424463780387327037251069E8L,
00074   Q3 = -2.019684072836541751428967854947019415698E7L,
00075   Q4 = 1.682912729190313538934190635536631941751E6L,
00076   Q5 = -9.615511549171441430850103489315371768998E4L,
00077   Q6 = 3.697714952261803935521187272204485251835E3L,
00078   Q7 = -8.802340681794263968892934703309274564037E1L,
00079   /* Q8 = 1.000000000000000000000000000000000000000E0 */
00080 /* C1 + C2 = ln 2 */
00081 
00082   C1 = 6.93145751953125E-1L,
00083   C2 = 1.428606820309417232121458176568075500134E-6L,
00084 /* ln (2^16384 * (1 - 2^-113)) */
00085   maxlog = 1.1356523406294143949491931077970764891253E4L,
00086 /* ln 2^-114 */
00087   minarg = -7.9018778583833765273564461846232128760607E1L, big = 2e307L;
00088 
00089 
00090 long double
00091 __expm1l (long double x)
00092 {
00093   long double px, qx, xx;
00094   int32_t ix, sign;
00095   ieee854_long_double_shape_type u;
00096   int k;
00097 
00098   /* Detect infinity and NaN.  */
00099   u.value = x;
00100   ix = u.parts32.w0;
00101   sign = ix & 0x80000000;
00102   ix &= 0x7fffffff;
00103   if (ix >= 0x7ff00000)
00104     {
00105       /* Infinity. */
00106       if (((ix & 0xfffff) | u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
00107        {
00108          if (sign)
00109            return -1.0L;
00110          else
00111            return x;
00112        }
00113       /* NaN. No invalid exception. */
00114       return x;
00115     }
00116 
00117   /* expm1(+- 0) = +- 0.  */
00118   if ((ix == 0) && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
00119     return x;
00120 
00121   /* Overflow.  */
00122   if (x > maxlog)
00123     return (big * big);
00124 
00125   /* Minimum value.  */
00126   if (x < minarg)
00127     return (4.0/big - 1.0L);
00128 
00129   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
00130   xx = C1 + C2;                    /* ln 2. */
00131   px = __floorl (0.5 + x / xx);
00132   k = px;
00133   /* remainder times ln 2 */
00134   x -= px * C1;
00135   x -= px * C2;
00136 
00137   /* Approximate exp(remainder ln 2).  */
00138   px = (((((((P7 * x
00139              + P6) * x
00140             + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
00141 
00142   qx = (((((((x
00143              + Q7) * x
00144             + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
00145 
00146   xx = x * x;
00147   qx = x + (0.5 * xx + xx * px / qx);
00148 
00149   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
00150 
00151   We have qx = exp(remainder ln 2) - 1, so
00152   exp(x) - 1 = 2^k (qx + 1) - 1
00153              = 2^k qx + 2^k - 1.  */
00154 
00155   px = __ldexpl (1.0L, k);
00156   x = px * qx + (px - 1.0);
00157   return x;
00158 }
00159 libm_hidden_def (__expm1l)
00160 long_double_symbol (libm, __expm1l, expm1l);