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