libquadmath/math/expm1q.c

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