BSD 4_3_Tahoe release
[unix-history] / usr / src / old / libom / log.c
CommitLineData
95f51977 1/* @(#)log.c 4.2 1/22/85 */
d571f5dd
SL
2
3/*
4 log returns the natural logarithm of its floating
5 point argument.
6
7 The coefficients are #2705 from Hart & Cheney. (19.38D)
8
9 It calls frexp.
10*/
11
12#include <errno.h>
13#include <math.h>
14
15int errno;
16double frexp();
17static double log2 = 0.693147180559945309e0;
18static double ln10 = 2.302585092994045684;
19static double sqrto2 = 0.707106781186547524e0;
20static double p0 = -.240139179559210510e2;
21static double p1 = 0.309572928215376501e2;
a52057f1 22static double p2 = -.963769093377840513e1;
d571f5dd
SL
23static double p3 = 0.421087371217979714e0;
24static double q0 = -.120069589779605255e2;
25static double q1 = 0.194809660700889731e2;
26static double q2 = -.891110902798312337e1;
27
28double
29log(arg)
30double arg;
31{
32 double x,z, zsq, temp;
33 int exp;
34
35 if(arg <= 0.) {
36 errno = EDOM;
37 return(-HUGE);
38 }
39 x = frexp(arg,&exp);
40 while(x<0.5) {
41 x = x*2;
42 exp = exp-1;
43 }
44 if(x<sqrto2) {
45 x = 2*x;
46 exp = exp-1;
47 }
48
49 z = (x-1)/(x+1);
50 zsq = z*z;
51
52 temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
53 temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0);
54 temp = temp*z + exp*log2;
55 return(temp);
56}
57
58double
59log10(arg)
60double arg;
61{
62
63 return(log(arg)/ln10);
64}