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
130
131
132
|
/* @(#)z_logarithm.c 1.0 98/08/13 */
/******************************************************************
* The following routines are coded directly from the algorithms
* and coefficients given in "Software Manual for the Elementary
* Functions" by William J. Cody, Jr. and William Waite, Prentice
* Hall, 1980.
******************************************************************/
/*
FUNCTION
<<log>>, <<logf>>, <<log10>>, <<log10f>>, <<logarithm>>, <<logarithmf>>---natural or base 10 logarithms
INDEX
log
INDEX
logf
INDEX
log10
INDEX
log10f
SYNOPSIS
#include <math.h>
double log(double <[x]>);
float logf(float <[x]>);
double log10(double <[x]>);
float log10f(float <[x]>);
DESCRIPTION
Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e
(where e is the base of the natural system of logarithms, 2.71828@dots{}) or
base 10.
<<log>> and <<logf>> are identical save for the return and argument types.
<<log10>> and <<log10f>> are identical save for the return and argument types.
RETURNS
Normally, returns the calculated value. When <[x]> is zero, the
returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>.
When <[x]> is negative, the returned value is <<-HUGE_VAL>> and
<<errno>> is set to <<EDOM>>.
PORTABILITY
<<log>> is ANSI. <<logf>> is an extension.
<<log10>> is ANSI. <<log10f>> is an extension.
*/
/******************************************************************
* Logarithm
*
* Input:
* x - floating point value
* ten - indicates base ten numbers
*
* Output:
* logarithm of x
*
* Description:
* This routine calculates logarithms.
*
*****************************************************************/
#include "fdlibm.h"
#include "zmath.h"
#ifndef _DOUBLE_IS_32BITS
static const double a[] = { -0.64124943423745581147e+02,
0.16383943563021534222e+02,
-0.78956112887481257267 };
static const double b[] = { -0.76949932108494879777e+03,
0.31203222091924532844e+03,
-0.35667977739034646171e+02 };
static const double C1 = 22713.0 / 32768.0;
static const double C2 = 1.428606820309417232e-06;
static const double C3 = 0.43429448190325182765;
double
logarithm (double x,
int ten)
{
int N;
double f, w, z;
/* Check for range and domain errors here. */
if (x == 0.0)
{
errno = ERANGE;
return (-z_infinity.d);
}
else if (x < 0.0)
{
errno = EDOM;
return (z_notanum.d);
}
else if (!isfinite(x))
{
if (isnan(x))
return (z_notanum.d);
else
return (z_infinity.d);
}
/* Get the exponent and mantissa where x = f * 2^N. */
f = frexp (x, &N);
z = f - 0.5;
if (f > __SQRT_HALF)
z = (z - 0.5) / (f * 0.5 + 0.5);
else
{
N--;
z /= (z * 0.5 + 0.5);
}
w = z * z;
/* Use Newton's method with 4 terms. */
z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]);
if (N != 0)
z = (N * C2 + z) + N * C1;
if (ten)
z *= C3;
return (z);
}
#endif /* _DOUBLE_IS_32BITS */
|
