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