64d5146139fd4acab98fb838e5871569d1c9812d
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-128 / s_log1pl.c
1 /*                                                      log1pl.c
2  *
3  *      Relative error logarithm
4  *      Natural logarithm of 1+x, 128-bit long double precision
5  *
6  *
7  *
8  * SYNOPSIS:
9  *
10  * long double x, y, log1pl();
11  *
12  * y = log1pl( x );
13  *
14  *
15  *
16  * DESCRIPTION:
17  *
18  * Returns the base e (2.718...) logarithm of 1+x.
19  *
20  * The argument 1+x is separated into its exponent and fractional
21  * parts.  If the exponent is between -1 and +1, the logarithm
22  * of the fraction is approximated by
23  *
24  *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
25  *
26  * Otherwise, setting  z = 2(w-1)/(w+1),
27  *
28  *     log(w) = z + z^3 P(z)/Q(z).
29  *
30  *
31  *
32  * ACCURACY:
33  *
34  *                      Relative error:
35  * arithmetic   domain     # trials      peak         rms
36  *    IEEE      -1, 8       100000      1.9e-34     4.3e-35
37  */
38
39 /* Copyright 2001 by Stephen L. Moshier 
40
41     This library is free software; you can redistribute it and/or
42     modify it under the terms of the GNU Lesser General Public
43     License as published by the Free Software Foundation; either
44     version 2.1 of the License, or (at your option) any later version.
45
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, write to the Free Software
53     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
54
55
56 #include "math.h"
57 #include "math_private.h"
58
59 /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
60  * 1/sqrt(2) <= 1+x < sqrt(2)
61  * Theoretical peak relative error = 5.3e-37,
62  * relative peak error spread = 2.3e-14
63  */
64 static const long double
65   P12 = 1.538612243596254322971797716843006400388E-6L,
66   P11 = 4.998469661968096229986658302195402690910E-1L,
67   P10 = 2.321125933898420063925789532045674660756E1L,
68   P9 = 4.114517881637811823002128927449878962058E2L,
69   P8 = 3.824952356185897735160588078446136783779E3L,
70   P7 = 2.128857716871515081352991964243375186031E4L,
71   P6 = 7.594356839258970405033155585486712125861E4L,
72   P5 = 1.797628303815655343403735250238293741397E5L,
73   P4 = 2.854829159639697837788887080758954924001E5L,
74   P3 = 3.007007295140399532324943111654767187848E5L,
75   P2 = 2.014652742082537582487669938141683759923E5L,
76   P1 = 7.771154681358524243729929227226708890930E4L,
77   P0 = 1.313572404063446165910279910527789794488E4L,
78   /* Q12 = 1.000000000000000000000000000000000000000E0L, */
79   Q11 = 4.839208193348159620282142911143429644326E1L,
80   Q10 = 9.104928120962988414618126155557301584078E2L,
81   Q9 = 9.147150349299596453976674231612674085381E3L,
82   Q8 = 5.605842085972455027590989944010492125825E4L,
83   Q7 = 2.248234257620569139969141618556349415120E5L,
84   Q6 = 6.132189329546557743179177159925690841200E5L,
85   Q5 = 1.158019977462989115839826904108208787040E6L,
86   Q4 = 1.514882452993549494932585972882995548426E6L,
87   Q3 = 1.347518538384329112529391120390701166528E6L,
88   Q2 = 7.777690340007566932935753241556479363645E5L,
89   Q1 = 2.626900195321832660448791748036714883242E5L,
90   Q0 = 3.940717212190338497730839731583397586124E4L;
91
92 /* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
93  * where z = 2(x-1)/(x+1)
94  * 1/sqrt(2) <= x < sqrt(2)
95  * Theoretical peak relative error = 1.1e-35,
96  * relative peak error spread 1.1e-9
97  */
98 static const long double
99   R5 = -8.828896441624934385266096344596648080902E-1L,
100   R4 = 8.057002716646055371965756206836056074715E1L,
101   R3 = -2.024301798136027039250415126250455056397E3L,
102   R2 = 2.048819892795278657810231591630928516206E4L,
103   R1 = -8.977257995689735303686582344659576526998E4L,
104   R0 = 1.418134209872192732479751274970992665513E5L,
105   /* S6 = 1.000000000000000000000000000000000000000E0L, */
106   S5 = -1.186359407982897997337150403816839480438E2L,
107   S4 = 3.998526750980007367835804959888064681098E3L,
108   S3 = -5.748542087379434595104154610899551484314E4L,
109   S2 = 4.001557694070773974936904547424676279307E5L,
110   S1 = -1.332535117259762928288745111081235577029E6L,
111   S0 = 1.701761051846631278975701529965589676574E6L;
112
113 /* C1 + C2 = ln 2 */
114 static const long double C1 = 6.93145751953125E-1L;
115 static const long double C2 = 1.428606820309417232121458176568075500134E-6L;
116
117 static const long double sqrth = 0.7071067811865475244008443621048490392848L;
118 /* ln (2^16384 * (1 - 2^-113)) */
119 static const long double maxlog = 1.1356523406294143949491931077970764891253E4L;
120 static const long double big = 2e4932L;
121 static const long double zero = 0.0L;
122
123 #if 1
124 /* Make sure these are prototyped.  */
125 long double frexpl (long double, int *);
126 long double ldexpl (long double, int);
127 #endif
128
129
130 long double
131 __log1pl (long double xm1)
132 {
133   long double x, y, z, r, s;
134   ieee854_long_double_shape_type u;
135   int32_t hx;
136   int e;
137
138   /* Test for NaN or infinity input. */
139   u.value = xm1;
140   hx = u.parts32.w0;
141   if (hx >= 0x7fff0000)
142     return xm1;
143
144   /* log1p(+- 0) = +- 0.  */
145   if (((hx & 0x7fffffff) == 0)
146       && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
147     return xm1;
148
149   x = xm1 + 1.0L;
150
151   /* log1p(-1) = -inf */
152   if (x <= 0.0L)
153     {
154       if (x == 0.0L)
155         return (-1.0L / (x - x));
156       else
157         return (zero / (x - x));
158     }
159
160   /* Separate mantissa from exponent.  */
161
162   /* Use frexp used so that denormal numbers will be handled properly.  */
163   x = frexpl (x, &e);
164
165   /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
166      where z = 2(x-1)/x+1).  */
167   if ((e > 2) || (e < -2))
168     {
169       if (x < sqrth)
170         {                       /* 2( 2x-1 )/( 2x+1 ) */
171           e -= 1;
172           z = x - 0.5L;
173           y = 0.5L * z + 0.5L;
174         }
175       else
176         {                       /*  2 (x-1)/(x+1)   */
177           z = x - 0.5L;
178           z -= 0.5L;
179           y = 0.5L * x + 0.5L;
180         }
181       x = z / y;
182       z = x * x;
183       r = ((((R5 * z
184               + R4) * z
185              + R3) * z
186             + R2) * z
187            + R1) * z
188         + R0;
189       s = (((((z
190                + S5) * z
191               + S4) * z
192              + S3) * z
193             + S2) * z
194            + S1) * z
195         + S0;
196       z = x * (z * r / s);
197       z = z + e * C2;
198       z = z + x;
199       z = z + e * C1;
200       return (z);
201     }
202
203
204   /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
205
206   if (x < sqrth)
207     {
208       e -= 1;
209       if (e != 0)
210         x = 2.0L * x - 1.0L;    /*  2x - 1  */
211       else
212         x = xm1;
213     }
214   else
215     {
216       if (e != 0)
217         x = x - 1.0L;
218       else
219         x = xm1;
220     }
221   z = x * x;
222   r = (((((((((((P12 * x
223                  + P11) * x
224                 + P10) * x
225                + P9) * x
226               + P8) * x
227              + P7) * x
228             + P6) * x
229            + P5) * x
230           + P4) * x
231          + P3) * x
232         + P2) * x
233        + P1) * x
234     + P0;
235   s = (((((((((((x
236                  + Q11) * x
237                 + Q10) * x
238                + Q9) * x
239               + Q8) * x
240              + Q7) * x
241             + Q6) * x
242            + Q5) * x
243           + Q4) * x
244          + Q3) * x
245         + Q2) * x
246        + Q1) * x
247     + Q0;
248   y = x * (z * r / s);
249   y = y + e * C2;
250   z = y - 0.5L * z;
251   z = z + x;
252   z = z + e * C1;
253   return (z);
254 }
255
256 weak_alias (__log1pl, log1pl)