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 <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 log1pl(long double xm1)
115 {
116 long double x, y, z;
117 int e;
118 
119 if( isnan(xm1) )
120 	return(xm1);
121 if( xm1 == INFINITY )
122 	return(xm1);
123 if(xm1 == 0.0)
124 	return(xm1);
125 
126 x = xm1 + 1.0L;
127 
128 /* Test for domain errors.  */
129 if( x <= 0.0L )
130 	{
131 	if( x == 0.0L )
132 		return( -INFINITY );
133 	else
134 		return( NAN );
135 	}
136 
137 /* Separate mantissa from exponent.
138    Use frexp so that denormal numbers will be handled properly.  */
139 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 if( (e > 2) || (e < -2) )
144 {
145 if( x < SQRTH )
146 	{ /* 2( 2x-1 )/( 2x+1 ) */
147 	e -= 1;
148 	z = x - 0.5L;
149 	y = 0.5L * z + 0.5L;
150 	}
151 else
152 	{ /*  2 (x-1)/(x+1)   */
153 	z = x - 0.5L;
154 	z -= 0.5L;
155 	y = 0.5L * x  + 0.5L;
156 	}
157 x = z / y;
158 z = x*x;
159 z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
160 z = z + e * C2;
161 z = z + x;
162 z = z + e * C1;
163 return( z );
164 }
165 
166 
167 /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
168 
169 if( x < SQRTH )
170 	{
171 	e -= 1;
172 	if (e != 0)
173 	  x = 2.0 * x - 1.0L;
174 	else
175 	  x = xm1;
176 	}
177 else
178 	{
179 	  if (e != 0)
180 	    x = x - 1.0L;
181 	  else
182 	    x = xm1;
183 	}
184 z = x*x;
185 y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
186 y = y + e * C2;
187 z = y - 0.5 * z;
188 z = z + x;
189 z = z + e * C1;
190 return( z );
191 }
192