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  }