1 /* $OpenBSD: s_log1pl.c,v 1.1 2011/07/06 00:02:42 martynas Exp $ */
2
3 /*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19 /* log1pl.c
20 *
21 * Relative error logarithm
22 * Natural logarithm of 1+x, 128-bit long double precision
23 *
24 *
25 *
26 * SYNOPSIS:
27 *
28 * long double x, y, log1pl();
29 *
30 * y = log1pl( x );
31 *
32 *
33 *
34 * DESCRIPTION:
35 *
36 * Returns the base e (2.718...) logarithm of 1+x.
37 *
38 * The argument 1+x is separated into its exponent and fractional
39 * parts. If the exponent is between -1 and +1, the logarithm
40 * of the fraction is approximated by
41 *
42 * log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
43 *
44 * Otherwise, setting z = 2(w-1)/(w+1),
45 *
46 * log(w) = z + z^3 P(z)/Q(z).
47 *
48 *
49 *
50 * ACCURACY:
51 *
52 * Relative error:
53 * arithmetic domain # trials peak rms
54 * IEEE -1, 8 100000 1.9e-34 4.3e-35
55 */
56
57
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 zero = 0.0L;
120
121 long double
log1pl(long double xm1)122 log1pl(long double xm1)
123 {
124 long double x, y, z, r, s;
125 ieee_quad_shape_type u;
126 int32_t hx;
127 int e;
128
129 /* Test for NaN or infinity input. */
130 u.value = xm1;
131 hx = u.parts32.mswhi;
132 if (hx >= 0x7fff0000)
133 return xm1 + xm1;
134
135 /* log1p(+- 0) = +- 0. */
136 if (((hx & 0x7fffffff) == 0)
137 && (u.parts32.mswlo | u.parts32.lswhi | u.parts32.lswlo) == 0)
138 return xm1;
139
140 x = xm1 + 1.0L;
141
142 /* log1p(-1) = -inf */
143 if (x <= 0.0L)
144 {
145 if (x == 0.0L)
146 return __math_divzerol(1);
147 else
148 return __math_invalidl(xm1);
149 }
150
151 /* Separate mantissa from exponent. */
152
153 /* Use frexp used so that denormal numbers will be handled properly. */
154 x = frexpl (x, &e);
155
156 /* Logarithm using log(x) = z + z^3 P(z^2)/Q(z^2),
157 where z = 2(x-1)/x+1). */
158 if ((e > 2) || (e < -2))
159 {
160 if (x < sqrth)
161 { /* 2( 2x-1 )/( 2x+1 ) */
162 e -= 1;
163 z = x - 0.5L;
164 y = 0.5L * z + 0.5L;
165 }
166 else
167 { /* 2 (x-1)/(x+1) */
168 z = x - 0.5L;
169 z -= 0.5L;
170 y = 0.5L * x + 0.5L;
171 }
172 x = z / y;
173 z = x * x;
174 r = ((((R5 * z
175 + R4) * z
176 + R3) * z
177 + R2) * z
178 + R1) * z
179 + R0;
180 s = (((((z
181 + S5) * z
182 + S4) * z
183 + S3) * z
184 + S2) * z
185 + S1) * z
186 + S0;
187 z = x * (z * r / s);
188 z = z + e * C2;
189 z = z + x;
190 z = z + e * C1;
191 return (z);
192 }
193
194
195 /* Logarithm using log(1+x) = x - .5x^2 + x^3 P(x)/Q(x). */
196
197 if (x < sqrth)
198 {
199 e -= 1;
200 if (e != 0)
201 x = 2.0L * x - 1.0L; /* 2x - 1 */
202 else
203 x = xm1;
204 }
205 else
206 {
207 if (e != 0)
208 x = x - 1.0L;
209 else
210 x = xm1;
211 }
212 z = x * x;
213 r = (((((((((((P12 * x
214 + P11) * x
215 + P10) * x
216 + P9) * x
217 + P8) * x
218 + P7) * x
219 + P6) * x
220 + P5) * x
221 + P4) * x
222 + P3) * x
223 + P2) * x
224 + P1) * x
225 + P0;
226 s = (((((((((((x
227 + Q11) * x
228 + Q10) * x
229 + Q9) * x
230 + Q8) * x
231 + Q7) * x
232 + Q6) * x
233 + Q5) * x
234 + Q4) * x
235 + Q3) * x
236 + Q2) * x
237 + Q1) * x
238 + Q0;
239 y = x * (z * r / s);
240 y = y + e * C2;
241 z = y - 0.5L * z;
242 z = z + x;
243 z = z + e * C1;
244 return (z);
245 }
246
