1*9936Ssam /* @(#)log.c 4.1 12/25/82 */ 2*9936Ssam 3*9936Ssam /* 4*9936Ssam log returns the natural logarithm of its floating 5*9936Ssam point argument. 6*9936Ssam 7*9936Ssam The coefficients are #2705 from Hart & Cheney. (19.38D) 8*9936Ssam 9*9936Ssam It calls frexp. 10*9936Ssam */ 11*9936Ssam 12*9936Ssam #include <errno.h> 13*9936Ssam #include <math.h> 14*9936Ssam 15*9936Ssam int errno; 16*9936Ssam double frexp(); 17*9936Ssam static double log2 = 0.693147180559945309e0; 18*9936Ssam static double ln10 = 2.302585092994045684; 19*9936Ssam static double sqrto2 = 0.707106781186547524e0; 20*9936Ssam static double p0 = -.240139179559210510e2; 21*9936Ssam static double p1 = 0.309572928215376501e2; 22*9936Ssam static double p2 = -.963769093368686593e1; 23*9936Ssam static double p3 = 0.421087371217979714e0; 24*9936Ssam static double q0 = -.120069589779605255e2; 25*9936Ssam static double q1 = 0.194809660700889731e2; 26*9936Ssam static double q2 = -.891110902798312337e1; 27*9936Ssam 28*9936Ssam double 29*9936Ssam log(arg) 30*9936Ssam double arg; 31*9936Ssam { 32*9936Ssam double x,z, zsq, temp; 33*9936Ssam int exp; 34*9936Ssam 35*9936Ssam if(arg <= 0.) { 36*9936Ssam errno = EDOM; 37*9936Ssam return(-HUGE); 38*9936Ssam } 39*9936Ssam x = frexp(arg,&exp); 40*9936Ssam while(x<0.5) { 41*9936Ssam x = x*2; 42*9936Ssam exp = exp-1; 43*9936Ssam } 44*9936Ssam if(x<sqrto2) { 45*9936Ssam x = 2*x; 46*9936Ssam exp = exp-1; 47*9936Ssam } 48*9936Ssam 49*9936Ssam z = (x-1)/(x+1); 50*9936Ssam zsq = z*z; 51*9936Ssam 52*9936Ssam temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0; 53*9936Ssam temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0); 54*9936Ssam temp = temp*z + exp*log2; 55*9936Ssam return(temp); 56*9936Ssam } 57*9936Ssam 58*9936Ssam double 59*9936Ssam log10(arg) 60*9936Ssam double arg; 61*9936Ssam { 62*9936Ssam 63*9936Ssam return(log(arg)/ln10); 64*9936Ssam } 65