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 <float.h>
  58 #include <math.h>
  59 #include <math_private.h>
  60
  61 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
  62    -.5 ln 2  <  x  <  .5 ln 2
  63    Theoretical peak relative error = 8.1e-36  */
  64
  65 static const long double
  66   P0 = 2.943520915569954073888921213330863757240E8L,
  67   P1 = -5.722847283900608941516165725053359168840E7L,
  68   P2 = 8.944630806357575461578107295909719817253E6L,
  69   P3 = -7.212432713558031519943281748462837065308E5L,
  70   P4 = 4.578962475841642634225390068461943438441E4L,
  71   P5 = -1.716772506388927649032068540558788106762E3L,
  72   P6 = 4.401308817383362136048032038528753151144E1L,
  73   P7 = -4.888737542888633647784737721812546636240E-1L,
  74   Q0 = 1.766112549341972444333352727998584753865E9L,
  75   Q1 = -7.848989743695296475743081255027098295771E8L,
  76   Q2 = 1.615869009634292424463780387327037251069E8L,
  77   Q3 = -2.019684072836541751428967854947019415698E7L,
  78   Q4 = 1.682912729190313538934190635536631941751E6L,
  79   Q5 = -9.615511549171441430850103489315371768998E4L,
  80   Q6 = 3.697714952261803935521187272204485251835E3L,
  81   Q7 = -8.802340681794263968892934703309274564037E1L,
  82   /* Q8 = 1.000000000000000000000000000000000000000E0 */
  83 /* C1 + C2 = ln 2 */
  84
  85   C1 = 6.93145751953125E-1L,
  86   C2 = 1.428606820309417232121458176568075500134E-6L,
  87 /* ln (2^16384 * (1 - 2^-113)) */
  88   maxlog = 1.1356523406294143949491931077970764891253E4L,
  89 /* ln 2^-114 */
  90   minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e4932L;
  91
  92
  93 long double
  94 __expm1l (long double x)
  95 {
  96   long double px, qx, xx;
  97   int32_t ix, sign;
  98   ieee854_long_double_shape_type u;
  99   int k;
 100
 101   /* Detect infinity and NaN.  */
 102   u.value = x;
 103   ix = u.parts32.w0;
 104   sign = ix & 0x80000000;
 105   ix &= 0x7fffffff;
 106   if (!sign && ix >= 0x40060000)
 107     {
 108       /* If num is positive and exp >= 6 use plain exp.  */
 109       return __expl (x);
 110     }
 111   if (ix >= 0x7fff0000)
 112     {
 113       /* Infinity. */
 114       if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 115         {
 116           if (sign)
 117             return -1.0L;
 118           else
 119             return x;
 120         }
 121       /* NaN. No invalid exception. */
 122       return x;
 123     }
 124
 125   /* expm1(+- 0) = +- 0.  */
 126   if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 127     return x;
 128
 129   /* Overflow.  */
 130   if (x > maxlog)
 131     {
 132       __set_errno (ERANGE);
 133       return (big * big);
 134     }
 135
 136   /* Minimum value.  */
 137   if (x < minarg)
 138     return (4.0/big - 1.0L);
 139
 140   /* Avoid internal underflow when result does not underflow, while
 141      ensuring underflow (without returning a zero of the wrong sign)
 142      when the result does underflow.  */
 143   if (fabsl (x) < 0x1p-113L)
 144     {
 145       if (fabsl (x) < LDBL_MIN)
 146         {
 147           long double force_underflow = x * x;
 148           math_force_eval (force_underflow);
 149         }
 150       return x;
 151     }
 152
 153   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
 154   xx = C1 + C2;                 /* ln 2. */
 155   px = __floorl (0.5 + x / xx);
 156   k = px;
 157   /* remainder times ln 2 */
 158   x -= px * C1;
 159   x -= px * C2;
 160
 161   /* Approximate exp(remainder ln 2).  */
 162   px = (((((((P7 * x
 163               + P6) * x
 164              + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
 165
 166   qx = (((((((x
 167               + Q7) * x
 168              + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
 169
 170   xx = x * x;
 171   qx = x + (0.5 * xx + xx * px / qx);
 172
 173   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 174
 175   We have qx = exp(remainder ln 2) - 1, so
 176   exp(x) - 1 = 2^k (qx + 1) - 1
 177              = 2^k qx + 2^k - 1.  */
 178
 179   px = __ldexpl (1.0L, k);
 180   x = px * qx + (px - 1.0);
 181   return x;
 182 }
 183 libm_hidden_def (__expm1l)
 184 weak_alias (__expm1l, expm1l)