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)25log(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)52log10(double arg) 53 { 54 55 if(arg <= 0) 56 return NaN(); 57 return log(arg) * ln10o1; 58 } 59