1 |
|
|
/* $OpenBSD: e_logl.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 |
|
|
/* logl.c |
20 |
|
|
* |
21 |
|
|
* Natural logarithm, long double precision |
22 |
|
|
* |
23 |
|
|
* |
24 |
|
|
* |
25 |
|
|
* SYNOPSIS: |
26 |
|
|
* |
27 |
|
|
* long double x, y, logl(); |
28 |
|
|
* |
29 |
|
|
* y = logl( x ); |
30 |
|
|
* |
31 |
|
|
* |
32 |
|
|
* |
33 |
|
|
* DESCRIPTION: |
34 |
|
|
* |
35 |
|
|
* Returns the base e (2.718...) logarithm of x. |
36 |
|
|
* |
37 |
|
|
* The argument is separated into its exponent and fractional |
38 |
|
|
* parts. If the exponent is between -1 and +1, the logarithm |
39 |
|
|
* of the fraction is approximated by |
40 |
|
|
* |
41 |
|
|
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x). |
42 |
|
|
* |
43 |
|
|
* Otherwise, setting z = 2(x-1)/x+1), |
44 |
|
|
* |
45 |
|
|
* log(x) = z + z**3 P(z)/Q(z). |
46 |
|
|
* |
47 |
|
|
* |
48 |
|
|
* |
49 |
|
|
* ACCURACY: |
50 |
|
|
* |
51 |
|
|
* Relative error: |
52 |
|
|
* arithmetic domain # trials peak rms |
53 |
|
|
* IEEE 0.5, 2.0 150000 8.71e-20 2.75e-20 |
54 |
|
|
* IEEE exp(+-10000) 100000 5.39e-20 2.34e-20 |
55 |
|
|
* |
56 |
|
|
* In the tests over the interval exp(+-10000), the logarithms |
57 |
|
|
* of the random arguments were uniformly distributed over |
58 |
|
|
* [-10000, +10000]. |
59 |
|
|
* |
60 |
|
|
* ERROR MESSAGES: |
61 |
|
|
* |
62 |
|
|
* log singularity: x = 0; returns -INFINITY |
63 |
|
|
* log domain: x < 0; returns NAN |
64 |
|
|
*/ |
65 |
|
|
|
66 |
|
|
#include <math.h> |
67 |
|
|
|
68 |
|
|
#include "math_private.h" |
69 |
|
|
|
70 |
|
|
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) |
71 |
|
|
* 1/sqrt(2) <= x < sqrt(2) |
72 |
|
|
* Theoretical peak relative error = 2.32e-20 |
73 |
|
|
*/ |
74 |
|
|
static long double P[] = { |
75 |
|
|
4.5270000862445199635215E-5L, |
76 |
|
|
4.9854102823193375972212E-1L, |
77 |
|
|
6.5787325942061044846969E0L, |
78 |
|
|
2.9911919328553073277375E1L, |
79 |
|
|
6.0949667980987787057556E1L, |
80 |
|
|
5.7112963590585538103336E1L, |
81 |
|
|
2.0039553499201281259648E1L, |
82 |
|
|
}; |
83 |
|
|
static long double Q[] = { |
84 |
|
|
/* 1.0000000000000000000000E0,*/ |
85 |
|
|
1.5062909083469192043167E1L, |
86 |
|
|
8.3047565967967209469434E1L, |
87 |
|
|
2.2176239823732856465394E2L, |
88 |
|
|
3.0909872225312059774938E2L, |
89 |
|
|
2.1642788614495947685003E2L, |
90 |
|
|
6.0118660497603843919306E1L, |
91 |
|
|
}; |
92 |
|
|
|
93 |
|
|
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2), |
94 |
|
|
* where z = 2(x-1)/(x+1) |
95 |
|
|
* 1/sqrt(2) <= x < sqrt(2) |
96 |
|
|
* Theoretical peak relative error = 6.16e-22 |
97 |
|
|
*/ |
98 |
|
|
|
99 |
|
|
static long double R[4] = { |
100 |
|
|
1.9757429581415468984296E-3L, |
101 |
|
|
-7.1990767473014147232598E-1L, |
102 |
|
|
1.0777257190312272158094E1L, |
103 |
|
|
-3.5717684488096787370998E1L, |
104 |
|
|
}; |
105 |
|
|
static long double S[4] = { |
106 |
|
|
/* 1.00000000000000000000E0L,*/ |
107 |
|
|
-2.6201045551331104417768E1L, |
108 |
|
|
1.9361891836232102174846E2L, |
109 |
|
|
-4.2861221385716144629696E2L, |
110 |
|
|
}; |
111 |
|
|
static const long double C1 = 6.9314575195312500000000E-1L; |
112 |
|
|
static const long double C2 = 1.4286068203094172321215E-6L; |
113 |
|
|
|
114 |
|
|
#define SQRTH 0.70710678118654752440L |
115 |
|
|
|
116 |
|
|
long double |
117 |
|
|
logl(long double x) |
118 |
|
|
{ |
119 |
|
|
long double y, z; |
120 |
|
398 |
int e; |
121 |
|
|
|
122 |
✗✓ |
199 |
if( isnan(x) ) |
123 |
|
|
return(x); |
124 |
✗✓ |
199 |
if( x == INFINITY ) |
125 |
|
|
return(x); |
126 |
|
|
/* Test for domain */ |
127 |
✓✓ |
199 |
if( x <= 0.0L ) |
128 |
|
|
{ |
129 |
✗✓ |
98 |
if( x == 0.0L ) |
130 |
|
|
return( -INFINITY ); |
131 |
|
|
else |
132 |
|
98 |
return( NAN ); |
133 |
|
|
} |
134 |
|
|
|
135 |
|
|
/* separate mantissa from exponent */ |
136 |
|
|
|
137 |
|
|
/* Note, frexp is used so that denormal numbers |
138 |
|
|
* will be handled properly. |
139 |
|
|
*/ |
140 |
|
101 |
x = frexpl( x, &e ); |
141 |
|
|
|
142 |
|
|
/* logarithm using log(x) = z + z**3 P(z)/Q(z), |
143 |
|
|
* where z = 2(x-1)/x+1) |
144 |
|
|
*/ |
145 |
✓✓ |
101 |
if( (e > 2) || (e < -2) ) |
146 |
|
|
{ |
147 |
✓✗ |
28 |
if( x < SQRTH ) |
148 |
|
|
{ /* 2( 2x-1 )/( 2x+1 ) */ |
149 |
|
28 |
e -= 1; |
150 |
|
28 |
z = x - 0.5L; |
151 |
|
28 |
y = 0.5L * z + 0.5L; |
152 |
|
28 |
} |
153 |
|
|
else |
154 |
|
|
{ /* 2 (x-1)/(x+1) */ |
155 |
|
|
z = x - 0.5L; |
156 |
|
|
z -= 0.5L; |
157 |
|
|
y = 0.5L * x + 0.5L; |
158 |
|
|
} |
159 |
|
28 |
x = z / y; |
160 |
|
28 |
z = x*x; |
161 |
|
28 |
z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) ); |
162 |
|
28 |
z = z + e * C2; |
163 |
|
28 |
z = z + x; |
164 |
|
28 |
z = z + e * C1; |
165 |
|
28 |
return( z ); |
166 |
|
|
} |
167 |
|
|
|
168 |
|
|
|
169 |
|
|
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */ |
170 |
|
|
|
171 |
✓✓ |
73 |
if( x < SQRTH ) |
172 |
|
|
{ |
173 |
|
45 |
e -= 1; |
174 |
|
45 |
x = ldexpl( x, 1 ) - 1.0L; /* 2x - 1 */ |
175 |
|
45 |
} |
176 |
|
|
else |
177 |
|
|
{ |
178 |
|
28 |
x = x - 1.0L; |
179 |
|
|
} |
180 |
|
73 |
z = x*x; |
181 |
|
73 |
y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) ); |
182 |
|
73 |
y = y + e * C2; |
183 |
|
73 |
z = y - ldexpl( z, -1 ); /* y - 0.5 * z */ |
184 |
|
|
/* Note, the sum of above terms does not exceed x/4, |
185 |
|
|
* so it contributes at most about 1/4 lsb to the error. |
186 |
|
|
*/ |
187 |
|
73 |
z = z + x; |
188 |
|
73 |
z = z + e * C1; /* This sum has an error of 1/2 lsb. */ |
189 |
|
73 |
return( z ); |
190 |
|
199 |
} |
191 |
|
|
DEF_STD(logl); |