Branch data Line data Source code
1 : : /* $OpenBSD: s_log1pl.c,v 1.3 2013/11/12 20:35:19 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, 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(x-1)/x+1),
45 : : *
46 : : * log(x) = 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.0, 9.0 100000 8.2e-20 2.5e-20
55 : : *
56 : : * ERROR MESSAGES:
57 : : *
58 : : * log singularity: x-1 = 0; returns -INFINITY
59 : : * log domain: x-1 < 0; returns NAN
60 : : */
61 : :
62 : : #include <openlibm_math.h>
63 : :
64 : : #include "math_private.h"
65 : :
66 : : /* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
67 : : * 1/sqrt(2) <= x < sqrt(2)
68 : : * Theoretical peak relative error = 2.32e-20
69 : : */
70 : :
71 : : static long double P[] = {
72 : : 4.5270000862445199635215E-5L,
73 : : 4.9854102823193375972212E-1L,
74 : : 6.5787325942061044846969E0L,
75 : : 2.9911919328553073277375E1L,
76 : : 6.0949667980987787057556E1L,
77 : : 5.7112963590585538103336E1L,
78 : : 2.0039553499201281259648E1L,
79 : : };
80 : : static long double Q[] = {
81 : : /* 1.0000000000000000000000E0,*/
82 : : 1.5062909083469192043167E1L,
83 : : 8.3047565967967209469434E1L,
84 : : 2.2176239823732856465394E2L,
85 : : 3.0909872225312059774938E2L,
86 : : 2.1642788614495947685003E2L,
87 : : 6.0118660497603843919306E1L,
88 : : };
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 = 6.16e-22
94 : : */
95 : :
96 : : static long double R[4] = {
97 : : 1.9757429581415468984296E-3L,
98 : : -7.1990767473014147232598E-1L,
99 : : 1.0777257190312272158094E1L,
100 : : -3.5717684488096787370998E1L,
101 : : };
102 : : static long double S[4] = {
103 : : /* 1.00000000000000000000E0L,*/
104 : : -2.6201045551331104417768E1L,
105 : : 1.9361891836232102174846E2L,
106 : : -4.2861221385716144629696E2L,
107 : : };
108 : : static const long double C1 = 6.9314575195312500000000E-1L;
109 : : static const long double C2 = 1.4286068203094172321215E-6L;
110 : :
111 : : #define SQRTH 0.70710678118654752440L
112 : :
113 : : long double
114 : 0 : log1pl(long double xm1)
115 : : {
116 : : long double x, y, z;
117 : : int e;
118 : :
119 [ # # ]: 0 : if( isnan(xm1) )
120 : 0 : return(xm1);
121 [ # # ]: 0 : if( xm1 == INFINITY )
122 : 0 : return(xm1);
123 [ # # ]: 0 : if(xm1 == 0.0)
124 : 0 : return(xm1);
125 : :
126 : 0 : x = xm1 + 1.0L;
127 : :
128 : : /* Test for domain errors. */
129 [ # # ]: 0 : if( x <= 0.0L )
130 : : {
131 [ # # ]: 0 : if( x == 0.0L )
132 : 0 : return( -INFINITY );
133 : : else
134 : 0 : return( NAN );
135 : : }
136 : :
137 : : /* Separate mantissa from exponent.
138 : : Use frexp so that denormal numbers will be handled properly. */
139 : 0 : x = frexpl( x, &e );
140 : :
141 : : /* logarithm using log(x) = z + z^3 P(z)/Q(z),
142 : : where z = 2(x-1)/x+1) */
143 [ # # # # ]: 0 : if( (e > 2) || (e < -2) )
144 : : {
145 [ # # ]: 0 : if( x < SQRTH )
146 : : { /* 2( 2x-1 )/( 2x+1 ) */
147 : 0 : e -= 1;
148 : 0 : z = x - 0.5L;
149 : 0 : y = 0.5L * z + 0.5L;
150 : : }
151 : : else
152 : : { /* 2 (x-1)/(x+1) */
153 : 0 : z = x - 0.5L;
154 : 0 : z -= 0.5L;
155 : 0 : y = 0.5L * x + 0.5L;
156 : : }
157 : 0 : x = z / y;
158 : 0 : z = x*x;
159 : 0 : z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
160 : 0 : z = z + e * C2;
161 : 0 : z = z + x;
162 : 0 : z = z + e * C1;
163 : 0 : return( z );
164 : : }
165 : :
166 : :
167 : : /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
168 : :
169 [ # # ]: 0 : if( x < SQRTH )
170 : : {
171 : 0 : e -= 1;
172 [ # # ]: 0 : if (e != 0)
173 : 0 : x = 2.0 * x - 1.0L;
174 : : else
175 : 0 : x = xm1;
176 : : }
177 : : else
178 : : {
179 [ # # ]: 0 : if (e != 0)
180 : 0 : x = x - 1.0L;
181 : : else
182 : 0 : x = xm1;
183 : : }
184 : 0 : z = x*x;
185 : 0 : y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
186 : 0 : y = y + e * C2;
187 : 0 : z = y - 0.5 * z;
188 : 0 : z = z + x;
189 : 0 : z = z + e * C1;
190 : 0 : return( z );
191 : : }
|
