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 <lib9.h> 11 12 #define log2 0.693147180559945309e0 13 #define ln10o1 .4342944819032518276511 14 #define sqrto2 0.707106781186547524e0 15 #define p0 -.240139179559210510e2 16 #define p1 0.309572928215376501e2 17 #define p2 -.963769093377840513e1 18 #define p3 0.421087371217979714e0 19 #define q0 -.120069589779605255e2 20 #define q1 0.194809660700889731e2 21 #define q2 -.891110902798312337e1 22 23 double log(double arg)24log(double arg) 25 { 26 double x, z, zsq, temp; 27 int exp; 28 29 if(arg <= 0) 30 return NaN(); 31 x = frexp(arg, &exp); 32 while(x < 0.5) { 33 x *= 2; 34 exp--; 35 } 36 if(x < sqrto2) { 37 x *= 2; 38 exp--; 39 } 40 41 z = (x-1) / (x+1); 42 zsq = z*z; 43 44 temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0; 45 temp = temp/(((zsq + q2)*zsq + q1)*zsq + q0); 46 temp = temp*z + exp*log2; 47 return temp; 48 } 49 50 double log10(double arg)51log10(double arg) 52 { 53 54 if(arg <= 0) 55 return NaN(); 56 return log(arg) * ln10o1; 57 } 58