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