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
129
|
/* log10f.c
*
* Common logarithm
*
*
*
* SYNOPSIS:
*
* float x, y, log10f();
*
* y = log10f( x );
*
*
*
* DESCRIPTION:
*
* Returns logarithm to the base 10 of x.
*
* The argument is separated into its exponent and fractional
* parts. 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 1.3e-7 3.4e-8
* IEEE 0, MAXNUMF 100000 1.3e-7 2.6e-8
*
* In the tests over the interval [0, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [-MAXL10, MAXL10].
*
* ERROR MESSAGES:
*
* log10f singularity: x = 0; returns -MAXL10
* log10f domain: x < 0; returns -MAXL10
* MAXL10 = 38.230809449325611792
*/
/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <math.h>
static char fname[] = {"log10"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
static float P[] = {
7.0376836292E-2,
-1.1514610310E-1,
1.1676998740E-1,
-1.2420140846E-1,
1.4249322787E-1,
-1.6668057665E-1,
2.0000714765E-1,
-2.4999993993E-1,
3.3333331174E-1
};
#define SQRTH 0.70710678118654752440
#define L102A 3.0078125E-1
#define L102B 2.48745663981195213739E-4
#define L10EA 4.3359375E-1
#define L10EB 7.00731903251827651129E-4
static float MAXL10 = 38.230809449325611792;
float frexpf(float, int *), polevlf(float, float *, int);
float log10f(float xx)
{
float x, y, z;
int e;
x = xx;
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
mtherr( fname, SING );
else
mtherr( fname, DOMAIN );
return( -MAXL10 );
}
/* separate mantissa from exponent */
x = frexpf( x, &e );
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x) */
if( x < SQRTH )
{
e -= 1;
x = 2.0*x - 1.0;
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
y = x * ( z * polevlf( x, P, 8 ) );
y = y - 0.5 * z; /* y - 0.5 * x**2 */
/* multiply log of fraction by log10(e)
* and base 2 exponent by log10(2)
*/
z = (x + y) * L10EB; /* accumulate terms in order of size */
z += y * L10EA;
z += x * L10EA;
x = e;
z += x * L102B;
z += x * L102A;
return( z );
}
|
