1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
/* logf.c
*
* Natural logarithm
*
*
*
* SYNOPSIS:
*
* float x, y, logf();
*
* y = logf( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 100000 7.6e-8 2.7e-8
* IEEE 1, MAXNUMF 100000 2.6e-8
*
* In the tests over the interval [1, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOGF].
*
* ERROR MESSAGES:
*
* logf singularity: x = 0; returns MINLOG
* logf domain: x < 0; returns MINLOG
*/
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Single precision natural logarithm
* test interval: [sqrt(2)/2, sqrt(2)]
* trials: 10000
* peak relative error: 7.1e-8
* rms relative error: 2.7e-8
*/
#include <math.h>
extern float MINLOGF, SQRTHF;
float frexpf( float, int * );
float logf( float xx )
{
register float y;
float x, z, fe;
int e;
x = xx;
fe = 0.0;
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
mtherr( "logf", SING );
else
mtherr( "logf", DOMAIN );
return( MINLOGF );
}
x = frexpf( x, &e );
if( x < SQRTHF )
{
e -= 1;
x = x + x - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
z = x * x;
/* 3.4e-9 */
/*
p = logfcof;
y = *p++ * x;
for( i=0; i<8; i++ )
{
y += *p++;
y *= x;
}
y *= z;
*/
y =
(((((((( 7.0376836292E-2 * x
- 1.1514610310E-1) * x
+ 1.1676998740E-1) * x
- 1.2420140846E-1) * x
+ 1.4249322787E-1) * x
- 1.6668057665E-1) * x
+ 2.0000714765E-1) * x
- 2.4999993993E-1) * x
+ 3.3333331174E-1) * x * z;
if( e )
{
fe = e;
y += -2.12194440e-4 * fe;
}
y += -0.5 * z; /* y - 0.5 x^2 */
z = x + y; /* ... + x */
if( e )
z += 0.693359375 * fe;
return( z );
}
|
