xref: /plan9/sys/src/libc/port/log.c (revision 3e12c5d1bb89fc02707907988834ef147769ddaf)
1 /*
2 	log returns the natural logarithm of its floating
3 	point argument.
4 
5 	The coefficients are #2705 from Hart & Cheney. (19.38D)
6 
7 	It calls frexp.
8 */
9 
10 #include <u.h>
11 #include <libc.h>
12 
13 #define	log2    0.693147180559945309e0
14 #define	ln10o1  .4342944819032518276511
15 #define	sqrto2  0.707106781186547524e0
16 #define	p0      -.240139179559210510e2
17 #define	p1      0.309572928215376501e2
18 #define	p2      -.963769093377840513e1
19 #define	p3      0.421087371217979714e0
20 #define	q0      -.120069589779605255e2
21 #define	q1      0.194809660700889731e2
22 #define	q2      -.891110902798312337e1
23 
24 double
log(double arg)25 log(double arg)
26 {
27 	double x, z, zsq, temp;
28 	int exp;
29 
30 	if(arg <= 0)
31 		return NaN();
32 	x = frexp(arg, &exp);
33 	while(x < 0.5) {
34 		x *= 2;
35 		exp--;
36 	}
37 	if(x < sqrto2) {
38 		x *= 2;
39 		exp--;
40 	}
41 
42 	z = (x-1) / (x+1);
43 	zsq = z*z;
44 
45 	temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
46 	temp = temp/(((zsq + q2)*zsq + q1)*zsq + q0);
47 	temp = temp*z + exp*log2;
48 	return temp;
49 }
50 
51 double
log10(double arg)52 log10(double arg)
53 {
54 
55 	if(arg <= 0)
56 		return NaN();
57 	return log(arg) * ln10o1;
58 }
59