GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libm/src/ld80/s_log1pl.c Lines: 0 41 0.0 %
Date: 2017-11-13 Branches: 0 20 0.0 %

Line Branch Exec Source
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);