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 Colombierlog(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 Colombierlog10(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