sysdeps/ieee754/ldbl-128ibm/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, write to the Free Software
  52     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  53
  54 #include <errno.h>
  55 #include "math.h"
  56 #include "math_private.h"
  57 #include <math_ldbl_opt.h>
  58
  59 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
  60    -.5 ln 2  <  x  <  .5 ln 2
  61    Theoretical peak relative error = 8.1e-36  */
  62
  63 static const long double
  64   P0 = 2.943520915569954073888921213330863757240E8L,
  65   P1 = -5.722847283900608941516165725053359168840E7L,
  66   P2 = 8.944630806357575461578107295909719817253E6L,
  67   P3 = -7.212432713558031519943281748462837065308E5L,
  68   P4 = 4.578962475841642634225390068461943438441E4L,
  69   P5 = -1.716772506388927649032068540558788106762E3L,
  70   P6 = 4.401308817383362136048032038528753151144E1L,
  71   P7 = -4.888737542888633647784737721812546636240E-1L,
  72   Q0 = 1.766112549341972444333352727998584753865E9L,
  73   Q1 = -7.848989743695296475743081255027098295771E8L,
  74   Q2 = 1.615869009634292424463780387327037251069E8L,
  75   Q3 = -2.019684072836541751428967854947019415698E7L,
  76   Q4 = 1.682912729190313538934190635536631941751E6L,
  77   Q5 = -9.615511549171441430850103489315371768998E4L,
  78   Q6 = 3.697714952261803935521187272204485251835E3L,
  79   Q7 = -8.802340681794263968892934703309274564037E1L,
  80   /* Q8 = 1.000000000000000000000000000000000000000E0 */
  81 /* C1 + C2 = ln 2 */
  82
  83   C1 = 6.93145751953125E-1L,
  84   C2 = 1.428606820309417232121458176568075500134E-6L,
  85 /* ln (2^16384 * (1 - 2^-113)) */
  86   maxlog = 1.1356523406294143949491931077970764891253E4L,
  87 /* ln 2^-114 */
  88   minarg = -7.9018778583833765273564461846232128760607E1L, big = 2e307L;
  89
  90
  91 long double
  92 __expm1l (long double x)
  93 {
  94   long double 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 (ix >= 0x7ff00000)
 105     {
 106       /* Infinity. */
 107       if (((ix & 0xfffff) | u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
 108         {
 109           if (sign)
 110             return -1.0L;
 111           else
 112             return x;
 113         }
 114       /* NaN. No invalid exception. */
 115       return x;
 116     }
 117
 118   /* expm1(+- 0) = +- 0.  */
 119   if ((ix == 0) && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
 120     return x;
 121
 122   /* Overflow.  */
 123   if (x > maxlog)
 124     {
 125       __set_errno (ERANGE);
 126       return (big * big);
 127     }
 128
 129   /* Minimum value.  */
 130   if (x < minarg)
 131     return (4.0/big - 1.0L);
 132
 133   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
 134   xx = C1 + C2;                 /* ln 2. */
 135   px = __floorl (0.5 + x / xx);
 136   k = px;
 137   /* remainder times ln 2 */
 138   x -= px * C1;
 139   x -= px * C2;
 140
 141   /* Approximate exp(remainder ln 2).  */
 142   px = (((((((P7 * x
 143               + P6) * x
 144              + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
 145
 146   qx = (((((((x
 147               + Q7) * x
 148              + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
 149
 150   xx = x * x;
 151   qx = x + (0.5 * xx + xx * px / qx);
 152
 153   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 154
 155   We have qx = exp(remainder ln 2) - 1, so
 156   exp(x) - 1 = 2^k (qx + 1) - 1
 157              = 2^k qx + 2^k - 1.  */
 158
 159   px = __ldexpl (1.0L, k);
 160   x = px * qx + (px - 1.0);
 161   return x;
 162 }
 163 libm_hidden_def (__expm1l)
 164 long_double_symbol (libm, __expm1l, expm1l);