libquadmath/math/expm1q.c

   1 /*                                                      expm1q.c
   2  *
   3  *      Exponential function, minus 1
   4  *      128-bit long double precision
   5  *
   6  *
   7  *
   8  * SYNOPSIS:
   9  *
  10  * long double x, y, expm1q();
  11  *
  12  * y = expm1q( 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 #include "quadmath-imp.h"
  55
  56 /* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
  57    -.5 ln 2  <  x  <  .5 ln 2
  58    Theoretical peak relative error = 8.1e-36  */
  59
  60 static const __float128
  61   P0 = 2.943520915569954073888921213330863757240E8Q,
  62   P1 = -5.722847283900608941516165725053359168840E7Q,
  63   P2 = 8.944630806357575461578107295909719817253E6Q,
  64   P3 = -7.212432713558031519943281748462837065308E5Q,
  65   P4 = 4.578962475841642634225390068461943438441E4Q,
  66   P5 = -1.716772506388927649032068540558788106762E3Q,
  67   P6 = 4.401308817383362136048032038528753151144E1Q,
  68   P7 = -4.888737542888633647784737721812546636240E-1Q,
  69   Q0 = 1.766112549341972444333352727998584753865E9Q,
  70   Q1 = -7.848989743695296475743081255027098295771E8Q,
  71   Q2 = 1.615869009634292424463780387327037251069E8Q,
  72   Q3 = -2.019684072836541751428967854947019415698E7Q,
  73   Q4 = 1.682912729190313538934190635536631941751E6Q,
  74   Q5 = -9.615511549171441430850103489315371768998E4Q,
  75   Q6 = 3.697714952261803935521187272204485251835E3Q,
  76   Q7 = -8.802340681794263968892934703309274564037E1Q,
  77   /* Q8 = 1.000000000000000000000000000000000000000E0 */
  78 /* C1 + C2 = ln 2 */
  79
  80   C1 = 6.93145751953125E-1Q,
  81   C2 = 1.428606820309417232121458176568075500134E-6Q,
  82 /* ln 2^-114 */
  83   minarg = -7.9018778583833765273564461846232128760607E1Q, big = 1e4932Q;
  84
  85
  86 __float128
  87 expm1q (__float128 x)
  88 {
  89   __float128 px, qx, xx;
  90   int32_t ix, sign;
  91   ieee854_float128 u;
  92   int k;
  93
  94   /* Detect infinity and NaN.  */
  95   u.value = x;
  96   ix = u.words32.w0;
  97   sign = ix & 0x80000000;
  98   ix &= 0x7fffffff;
  99   if (!sign && ix >= 0x40060000)
 100     {
 101       /* If num is positive and exp >= 6 use plain exp.  */
 102       return expq (x);
 103     }
 104   if (ix >= 0x7fff0000)
 105     {
 106       /* Infinity (which must be negative infinity). */
 107       if (((ix & 0xffff) | u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
 108         return -1;
 109       /* NaN.  Invalid exception if signaling.  */
 110       return x + x;
 111     }
 112
 113   /* expm1(+- 0) = +- 0.  */
 114   if ((ix == 0) && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
 115     return x;
 116
 117   /* Minimum value.  */
 118   if (x < minarg)
 119     return (4.0/big - 1);
 120
 121   /* Avoid internal underflow when result does not underflow, while
 122      ensuring underflow (without returning a zero of the wrong sign)
 123      when the result does underflow.  */
 124   if (fabsq (x) < 0x1p-113Q)
 125     {
 126       math_check_force_underflow (x);
 127       return x;
 128     }
 129
 130   /* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
 131   xx = C1 + C2;                 /* ln 2. */
 132   px = floorq (0.5 + x / xx);
 133   k = px;
 134   /* remainder times ln 2 */
 135   x -= px * C1;
 136   x -= px * C2;
 137
 138   /* Approximate exp(remainder ln 2).  */
 139   px = (((((((P7 * x
 140               + P6) * x
 141              + P5) * x + P4) * x + P3) * x + P2) * x + P1) * x + P0) * x;
 142
 143   qx = (((((((x
 144               + Q7) * x
 145              + Q6) * x + Q5) * x + Q4) * x + Q3) * x + Q2) * x + Q1) * x + Q0;
 146
 147   xx = x * x;
 148   qx = x + (0.5 * xx + xx * px / qx);
 149
 150   /* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
 151
 152   We have qx = exp(remainder ln 2) - 1, so
 153   exp(x) - 1 = 2^k (qx + 1) - 1
 154              = 2^k qx + 2^k - 1.  */
 155
 156   px = ldexpq (1, k);
 157   x = px * qx + (px - 1.0);
 158   return x;
 159 }