1 |
|
|
/* $OpenBSD: s_log1pl.c,v 1.5 2017/01/21 08:29:13 krw 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 |
|
|
DEF_STD(log1pl); |