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 <math.h> 11 #include <errno.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 errno = (arg==0)? ERANGE : EDOM; 32 return -HUGE_VAL; 33 } 34 x = frexp(arg, &exp); 35 while(x < 0.5) { 36 x *= 2; 37 exp--; 38 } 39 if(x < sqrto2) { 40 x *= 2; 41 exp--; 42 } 43 44 z = (x-1) / (x+1); 45 zsq = z*z; 46 47 temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0; 48 temp = temp/(((zsq + q2)*zsq + q1)*zsq + q0); 49 temp = temp*z + exp*log2; 50 return temp; 51 } 52 53 double log10(double arg)54log10(double arg) 55 { 56 57 if(arg <= 0) { 58 errno = (arg==0)? ERANGE : EDOM; 59 return -HUGE_VAL; 60 } 61 return log(arg) * ln10o1; 62 } 63