1 
2 /* @(#)z_logarithm.c 1.0 98/08/13 */
3 /******************************************************************
4  * The following routines are coded directly from the algorithms
5  * and coefficients given in "Software Manual for the Elementary
6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
7  * Hall, 1980.
8  ******************************************************************/
9 
10 /*
11 FUNCTION
12        <<log>>, <<logf>>, <<log10>>, <<log10f>>, <<logarithm>>, <<logarithmf>>---natural or base 10 logarithms
13 
14 INDEX
15     log
16 INDEX
17     logf
18 INDEX
19     log10
20 INDEX
21     log10f
22 
23 ANSI_SYNOPSIS
24        #include <math.h>
25        double log(double <[x]>);
26        float logf(float <[x]>);
27        double log10(double <[x]>);
28        float log10f(float <[x]>);
29 
30 TRAD_SYNOPSIS
31        #include <math.h>
32        double log(<[x]>);
33        double <[x]>;
34 
35        float logf(<[x]>);
36        float <[x]>;
37 
38        double log10(<[x]>);
39        double <[x]>;
40 
41        float log10f(<[x]>);
42        float <[x]>;
43 
44 DESCRIPTION
45 Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e
46 (where e is the base of the natural system of logarithms, 2.71828@dots{}) or
47 base 10.
48 <<log>> and <<logf>> are identical save for the return and argument types.
49 <<log10>> and <<log10f>> are identical save for the return and argument types.
50 
51 RETURNS
52 Normally, returns the calculated value.  When <[x]> is zero, the
53 returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>.
54 When <[x]> is negative, the returned value is <<-HUGE_VAL>> and
55 <<errno>> is set to <<EDOM>>.  You can control the error behavior via
56 <<matherr>>.
57 
58 PORTABILITY
59 <<log>> is ANSI. <<logf>> is an extension.
60 
61 <<log10>> is ANSI. <<log10f>> is an extension.
62 */
63 
64 
65 /******************************************************************
66  * Logarithm
67  *
68  * Input:
69  *   x - floating point value
70  *   ten - indicates base ten numbers
71  *
72  * Output:
73  *   logarithm of x
74  *
75  * Description:
76  *   This routine calculates logarithms.
77  *
78  *****************************************************************/
79 
80 #include "fdlibm.h"
81 #include "zmath.h"
82 
83 #ifndef _DOUBLE_IS_32BITS
84 
85 static const double a[] = { -0.64124943423745581147e+02,
86                              0.16383943563021534222e+02,
87                             -0.78956112887481257267 };
88 static const double b[] = { -0.76949932108494879777e+03,
89                              0.31203222091924532844e+03,
90                             -0.35667977739034646171e+02 };
91 static const double C1 =  22713.0 / 32768.0;
92 static const double C2 =  1.428606820309417232e-06;
93 static const double C3 =  0.43429448190325182765;
94 
95 double
96 _DEFUN (logarithm, (double, int),
97         double x _AND
98         int ten)
99 {
100   int N;
101   double f, w, z;
102 
103   /* Check for range and domain errors here. */
104   if (x == 0.0)
105     {
106       errno = ERANGE;
107       return (-z_infinity.d);
108     }
109   else if (x < 0.0)
110     {
111       errno = EDOM;
112       return (z_notanum.d);
113     }
114   else if (!isfinite(x))
115     {
116       if (isnan(x))
117         return (z_notanum.d);
118       else
119         return (z_infinity.d);
120     }
121 
122   /* Get the exponent and mantissa where x = f * 2^N. */
123   f = frexp (x, &N);
124 
125   z = f - 0.5;
126 
127   if (f > __SQRT_HALF)
128     z = (z - 0.5) / (f * 0.5 + 0.5);
129   else
130     {
131       N--;
132       z /= (z * 0.5 + 0.5);
133     }
134   w = z * z;
135 
136   /* Use Newton's method with 4 terms. */
137   z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]);
138 
139   if (N != 0)
140     z = (N * C2 + z) + N * C1;
141 
142   if (ten)
143     z *= C3;
144 
145   return (z);
146 }
147 
148 #endif /* _DOUBLE_IS_32BITS */
149