184447Sdfr/*							log1pq.c
284447Sdfr *
384447Sdfr *      Relative error logarithm
484447Sdfr *	Natural logarithm of 1+x, 128-bit long double precision
584447Sdfr *
684447Sdfr *
784447Sdfr *
884447Sdfr * SYNOPSIS:
984447Sdfr *
1084447Sdfr * long double x, y, log1pq();
1184447Sdfr *
1284447Sdfr * y = log1pq( x );
1384447Sdfr *
1484447Sdfr *
1584447Sdfr *
1684447Sdfr * DESCRIPTION:
1784447Sdfr *
1884447Sdfr * Returns the base e (2.718...) logarithm of 1+x.
1984447Sdfr *
2084447Sdfr * The argument 1+x is separated into its exponent and fractional
2184447Sdfr * parts.  If the exponent is between -1 and +1, the logarithm
2284447Sdfr * of the fraction is approximated by
2384447Sdfr *
2484447Sdfr *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
2584447Sdfr *
2684447Sdfr * Otherwise, setting  z = 2(w-1)/(w+1),
2784447Sdfr *
2884447Sdfr *     log(w) = z + z^3 P(z)/Q(z).
2984447Sdfr *
3084447Sdfr *
3184447Sdfr *
32193530Sjkim * ACCURACY:
33193530Sjkim *
3484447Sdfr *                      Relative error:
3584447Sdfr * arithmetic   domain     # trials      peak         rms
3684447Sdfr *    IEEE      -1, 8       100000      1.9e-34     4.3e-35
3784447Sdfr */
3884447Sdfr
3984447Sdfr/* Copyright 2001 by Stephen L. Moshier
4084447Sdfr
4184447Sdfr    This library is free software; you can redistribute it and/or
4284447Sdfr    modify it under the terms of the GNU Lesser General Public
4384447Sdfr    License as published by the Free Software Foundation; either
4484447Sdfr    version 2.1 of the License, or (at your option) any later version.
4584447Sdfr
46    This library is distributed in the hope that it will be useful,
47    but WITHOUT ANY WARRANTY; without even the implied warranty of
48    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
49    Lesser General Public License for more details.
50
51    You should have received a copy of the GNU Lesser General Public
52    License along with this library; if not, see
53    <http://www.gnu.org/licenses/>.  */
54
55#include "quadmath-imp.h"
56
57/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
58 * 1/sqrt(2) <= 1+x < sqrt(2)
59 * Theoretical peak relative error = 5.3e-37,
60 * relative peak error spread = 2.3e-14
61 */
62static const __float128
63  P12 = 1.538612243596254322971797716843006400388E-6Q,
64  P11 = 4.998469661968096229986658302195402690910E-1Q,
65  P10 = 2.321125933898420063925789532045674660756E1Q,
66  P9 = 4.114517881637811823002128927449878962058E2Q,
67  P8 = 3.824952356185897735160588078446136783779E3Q,
68  P7 = 2.128857716871515081352991964243375186031E4Q,
69  P6 = 7.594356839258970405033155585486712125861E4Q,
70  P5 = 1.797628303815655343403735250238293741397E5Q,
71  P4 = 2.854829159639697837788887080758954924001E5Q,
72  P3 = 3.007007295140399532324943111654767187848E5Q,
73  P2 = 2.014652742082537582487669938141683759923E5Q,
74  P1 = 7.771154681358524243729929227226708890930E4Q,
75  P0 = 1.313572404063446165910279910527789794488E4Q,
76  /* Q12 = 1.000000000000000000000000000000000000000E0L, */
77  Q11 = 4.839208193348159620282142911143429644326E1Q,
78  Q10 = 9.104928120962988414618126155557301584078E2Q,
79  Q9 = 9.147150349299596453976674231612674085381E3Q,
80  Q8 = 5.605842085972455027590989944010492125825E4Q,
81  Q7 = 2.248234257620569139969141618556349415120E5Q,
82  Q6 = 6.132189329546557743179177159925690841200E5Q,
83  Q5 = 1.158019977462989115839826904108208787040E6Q,
84  Q4 = 1.514882452993549494932585972882995548426E6Q,
85  Q3 = 1.347518538384329112529391120390701166528E6Q,
86  Q2 = 7.777690340007566932935753241556479363645E5Q,
87  Q1 = 2.626900195321832660448791748036714883242E5Q,
88  Q0 = 3.940717212190338497730839731583397586124E4Q;
89
90/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
91 * where z = 2(x-1)/(x+1)
92 * 1/sqrt(2) <= x < sqrt(2)
93 * Theoretical peak relative error = 1.1e-35,
94 * relative peak error spread 1.1e-9
95 */
96static const __float128
97  R5 = -8.828896441624934385266096344596648080902E-1Q,
98  R4 = 8.057002716646055371965756206836056074715E1Q,
99  R3 = -2.024301798136027039250415126250455056397E3Q,
100  R2 = 2.048819892795278657810231591630928516206E4Q,
101  R1 = -8.977257995689735303686582344659576526998E4Q,
102  R0 = 1.418134209872192732479751274970992665513E5Q,
103  /* S6 = 1.000000000000000000000000000000000000000E0L, */
104  S5 = -1.186359407982897997337150403816839480438E2Q,
105  S4 = 3.998526750980007367835804959888064681098E3Q,
106  S3 = -5.748542087379434595104154610899551484314E4Q,
107  S2 = 4.001557694070773974936904547424676279307E5Q,
108  S1 = -1.332535117259762928288745111081235577029E6Q,
109  S0 = 1.701761051846631278975701529965589676574E6Q;
110
111/* C1 + C2 = ln 2 */
112static const __float128 C1 = 6.93145751953125E-1Q;
113static const __float128 C2 = 1.428606820309417232121458176568075500134E-6Q;
114
115static const __float128 sqrth = 0.7071067811865475244008443621048490392848Q;
116/* ln (2^16384 * (1 - 2^-113)) */
117static const __float128 zero = 0;
118
119__float128
120log1pq (__float128 xm1)
121{
122  __float128 x, y, z, r, s;
123  ieee854_float128 u;
124  int32_t hx;
125  int e;
126
127  /* Test for NaN or infinity input. */
128  u.value = xm1;
129  hx = u.words32.w0;
130  if ((hx & 0x7fffffff) >= 0x7fff0000)
131    return xm1 + fabsq (xm1);
132
133  /* log1p(+- 0) = +- 0.  */
134  if (((hx & 0x7fffffff) == 0)
135      && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
136    return xm1;
137
138  if ((hx & 0x7fffffff) < 0x3f8e0000)
139    {
140      math_check_force_underflow (xm1);
141      if ((int) xm1 == 0)
142	return xm1;
143    }
144
145  if (xm1 >= 0x1p113Q)
146    x = xm1;
147  else
148    x = xm1 + 1;
149
150  /* log1p(-1) = -inf */
151  if (x <= 0)
152    {
153      if (x == 0)
154	return (-1 / zero);  /* log1p(-1) = -inf */
155      else
156	return (zero / (x - x));
157    }
158
159  /* Separate mantissa from exponent.  */
160
161  /* Use frexp used so that denormal numbers will be handled properly.  */
162  x = frexpq (x, &e);
163
164  /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
165     where z = 2(x-1)/x+1).  */
166  if ((e > 2) || (e < -2))
167    {
168      if (x < sqrth)
169	{			/* 2( 2x-1 )/( 2x+1 ) */
170	  e -= 1;
171	  z = x - 0.5Q;
172	  y = 0.5Q * z + 0.5Q;
173	}
174      else
175	{			/*  2 (x-1)/(x+1)   */
176	  z = x - 0.5Q;
177	  z -= 0.5Q;
178	  y = 0.5Q * x + 0.5Q;
179	}
180      x = z / y;
181      z = x * x;
182      r = ((((R5 * z
183	      + R4) * z
184	     + R3) * z
185	    + R2) * z
186	   + R1) * z
187	+ R0;
188      s = (((((z
189	       + S5) * z
190	      + S4) * z
191	     + S3) * z
192	    + S2) * z
193	   + S1) * z
194	+ S0;
195      z = x * (z * r / s);
196      z = z + e * C2;
197      z = z + x;
198      z = z + e * C1;
199      return (z);
200    }
201
202
203  /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
204
205  if (x < sqrth)
206    {
207      e -= 1;
208      if (e != 0)
209	x = 2 * x - 1;	/*  2x - 1  */
210      else
211	x = xm1;
212    }
213  else
214    {
215      if (e != 0)
216	x = x - 1;
217      else
218	x = xm1;
219    }
220  z = x * x;
221  r = (((((((((((P12 * x
222		 + P11) * x
223		+ P10) * x
224	       + P9) * x
225	      + P8) * x
226	     + P7) * x
227	    + P6) * x
228	   + P5) * x
229	  + P4) * x
230	 + P3) * x
231	+ P2) * x
232       + P1) * x
233    + P0;
234  s = (((((((((((x
235		 + Q11) * x
236		+ Q10) * x
237	       + Q9) * x
238	      + Q8) * x
239	     + Q7) * x
240	    + Q6) * x
241	   + Q5) * x
242	  + Q4) * x
243	 + Q3) * x
244	+ Q2) * x
245       + Q1) * x
246    + Q0;
247  y = x * (z * r / s);
248  y = y + e * C2;
249  z = y - 0.5Q * z;
250  z = z + x;
251  z = z + e * C1;
252  return (z);
253}
254