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 _Float128
  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^-114 */
  88   minarg = -7.9018778583833765273564461846232128760607E1L, big = 1e4932L;
  89
  90
  91 _Float128
  92 __expm1l (_Float128 x)
  93 {
  94   _Float128 px, qx, xx;
  95   int32_t ix, sign;
  96   ieee854_long_double_shape_type u;
  97   int k;
  98
  99   /* Detect infinity and NaN.  */
 100   u.value = x;
 101   ix = u.parts32.w0;
 102   sign = ix & 0x80000000;
 103   ix &= 0x7fffffff;
 104   if (!sign && ix >= 0x40060000)
 105     {
 106       /* If num is positive and exp >= 6 use plain exp.  */
 107       return __expl (x);
 108     }
 109   if (ix >= 0x7fff0000)
 110     {
 111       /* Infinity (which must be negative infinity). */
 112       if (((ix & 0xffff) | u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 113         return -1.0L;
 114       /* NaN.  Invalid exception if signaling.  */
 115       return x + x;
 116     }
 117
 118   /* expm1(+- 0) = +- 0.  */
 119   if ((ix == 0) && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 120     return x;
 121
 122   /* Minimum value.  */
 123   if (x < minarg)
 124     return (4.0/big - 1.0L);
 125
 126   /* Avoid internal underflow when result does not underflow, while
 127      ensuring underflow (without returning a zero of the wrong sign)
 128      when the result does underflow.  */
 129   if (fabsl (x) < 0x1p-113L)
 130     {
 131       math_check_force_underflow (x);
 132       return x;
 133     }
 134
 135   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
 136   xx = C1 + C2;                 /* ln 2. */
 137   px = __floorl (0.5 + x / xx);
 138   k = px;
 139   /* remainder times ln 2 */
 140   x -= px * C1;
 141   x -= px * C2;
 142
 143   /* Approximate exp(remainder ln 2).  */
 144   px = (((((((P7 * x
 145               + P6) * x
 146              + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
 147
 148   qx = (((((((x
 149               + Q7) * x
 150              + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
 151
 152   xx = x * x;
 153   qx = x + (0.5 * xx + xx * px / qx);
 154
 155   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 156
 157   We have qx = exp(remainder ln 2) - 1, so
 158   exp(x) - 1 = 2^k (qx + 1) - 1
 159              = 2^k qx + 2^k - 1.  */
 160
 161   px = __ldexpl (1.0L, k);
 162   x = px * qx + (px - 1.0);
 163   return x;
 164 }
 165 libm_hidden_def (__expm1l)
 166 weak_alias (__expm1l, expm1l)