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