sysdeps/ieee754/ldbl-128/s_expm1l.c

   1 /*                                                      expm1l.c
   2  *
   3  *      Exponential function, minus 1
   4  *      128-bit long double precision
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, expm1l();
  11  *
  12  * y = expm1l( x );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Returns e (2.71828...) raised to the x power, minus one.
  19  *
  20  * Range reduction is accomplished by separating the argument
  21  * into an integer k and fraction f such that
  22  *
  23  *     x    k  f
  24  *    e  = 2  e.
  25  *
  26  * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
  27  * in the basic range [-0.5 ln 2, 0.5 ln 2].
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain     # trials      peak         rms
  34  *    IEEE    -79,+MAXLOG    100,000     1.7e-34     4.5e-35
  35  *
  36  */
  37
  38 /* Copyright 2001 by Stephen L. Moshier
  39
  40     This library is free software; you can redistribute it and/or
  41     modify it under the terms of the GNU Lesser General Public
  42     License as published by the Free Software Foundation; either
  43     version 2.1 of the License, or (at your option) any later version.
  44
  45     This library is distributed in the hope that it will be useful,
  46     but WITHOUT ANY WARRANTY; without even the implied warranty of
  47     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  48     Lesser General Public License for more details.
  49
  50     You should have received a copy of the GNU Lesser General Public
  51     License along with this library; if not, see
  52     <http://www.gnu.org/licenses/>.  */
  53
  54
  55
  56 #include <errno.h>
  57 #include <math.h>
  58 #include <math_private.h>
  59
  60 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
  61    -.5 ln 2  <  x  <  .5 ln 2
  62    Theoretical peak relative error = 8.1e-36  */
  63
  64 static const long double
  65   P0 = 2.943520915569954073888921213330863757240E8L,
  66   P1 = -5.722847283900608941516165725053359168840E7L,
  67   P2 = 8.944630806357575461578107295909719817253E6L,
  68   P3 = -7.212432713558031519943281748462837065308E5L,
  69   P4 = 4.578962475841642634225390068461943438441E4L,
  70   P5 = -1.716772506388927649032068540558788106762E3L,
  71   P6 = 4.401308817383362136048032038528753151144E1L,
  72   P7 = -4.888737542888633647784737721812546636240E-1L,
  73   Q0 = 1.766112549341972444333352727998584753865E9L,
  74   Q1 = -7.848989743695296475743081255027098295771E8L,
  75   Q2 = 1.615869009634292424463780387327037251069E8L,
  76   Q3 = -2.019684072836541751428967854947019415698E7L,
  77   Q4 = 1.682912729190313538934190635536631941751E6L,
  78   Q5 = -9.615511549171441430850103489315371768998E4L,
  79   Q6 = 3.697714952261803935521187272204485251835E3L,
  80   Q7 = -8.802340681794263968892934703309274564037E1L,
  81   /* Q8 = 1.000000000000000000000000000000000000000E0 */
  82 /* C1 + C2 = ln 2 */
  83
  84   C1 = 6.93145751953125E-1L,
  85   C2 = 1.428606820309417232121458176568075500134E-6L,
  86 /* ln (2^16384 * (1 - 2^-113)) */
  87   maxlog = 1.1356523406294143949491931077970764891253E4L,
  88 /* ln 2^-114 */
  89   minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e4932L;
  90
  91
  92 long double
  93 __expm1l (long double x)
  94 {
  95   long double px, qx, xx;
  96   int32_t ix, sign;
  97   ieee854_long_double_shape_type u;
  98   int k;
  99
 100   /* Detect infinity and NaN.  */
 101   u.value = x;
 102   ix = u.parts32.w0;
 103   sign = ix & 0x80000000;
 104   ix &= 0x7fffffff;
 105   if (!sign && ix >= 0x40060000)
 106     {
 107       /* If num is positive and exp >= 6 use plain exp.  */
 108       return __expl (x);
 109     }
 110   if (ix >= 0x7fff0000)
 111     {
 112       /* Infinity. */
 113       if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 114         {
 115           if (sign)
 116             return -1.0L;
 117           else
 118             return x;
 119         }
 120       /* NaN. No invalid exception. */
 121       return x;
 122     }
 123
 124   /* expm1(+- 0) = +- 0.  */
 125   if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 126     return x;
 127
 128   /* Overflow.  */
 129   if (x > maxlog)
 130     {
 131       __set_errno (ERANGE);
 132       return (big * big);
 133     }
 134
 135   /* Minimum value.  */
 136   if (x < minarg)
 137     return (4.0/big - 1.0L);
 138
 139   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
 140   xx = C1 + C2;                 /* ln 2. */
 141   px = __floorl (0.5 + x / xx);
 142   k = px;
 143   /* remainder times ln 2 */
 144   x -= px * C1;
 145   x -= px * C2;
 146
 147   /* Approximate exp(remainder ln 2).  */
 148   px = (((((((P7 * x
 149               + P6) * x
 150              + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
 151
 152   qx = (((((((x
 153               + Q7) * x
 154              + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
 155
 156   xx = x * x;
 157   qx = x + (0.5 * xx + xx * px / qx);
 158
 159   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 160
 161   We have qx = exp(remainder ln 2) - 1, so
 162   exp(x) - 1 = 2^k (qx + 1) - 1
 163              = 2^k qx + 2^k - 1.  */
 164
 165   px = __ldexpl (1.0L, k);
 166   x = px * qx + (px - 1.0);
 167   return x;
 168 }
 169 libm_hidden_def (__expm1l)
 170 weak_alias (__expm1l, expm1l)