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
60static 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
87expm1q (__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}
160