1*17851Sralph /* @(#)log.c 4.2 01/22/85 */ 29936Ssam 39936Ssam /* 49936Ssam log returns the natural logarithm of its floating 59936Ssam point argument. 69936Ssam 79936Ssam The coefficients are #2705 from Hart & Cheney. (19.38D) 89936Ssam 99936Ssam It calls frexp. 109936Ssam */ 119936Ssam 129936Ssam #include <errno.h> 139936Ssam #include <math.h> 149936Ssam 159936Ssam int errno; 169936Ssam double frexp(); 179936Ssam static double log2 = 0.693147180559945309e0; 189936Ssam static double ln10 = 2.302585092994045684; 199936Ssam static double sqrto2 = 0.707106781186547524e0; 209936Ssam static double p0 = -.240139179559210510e2; 219936Ssam static double p1 = 0.309572928215376501e2; 22*17851Sralph static double p2 = -.963769093377840513e1; 239936Ssam static double p3 = 0.421087371217979714e0; 249936Ssam static double q0 = -.120069589779605255e2; 259936Ssam static double q1 = 0.194809660700889731e2; 269936Ssam static double q2 = -.891110902798312337e1; 279936Ssam 289936Ssam double log(arg)299936Ssamlog(arg) 309936Ssam double arg; 319936Ssam { 329936Ssam double x,z, zsq, temp; 339936Ssam int exp; 349936Ssam 359936Ssam if(arg <= 0.) { 369936Ssam errno = EDOM; 379936Ssam return(-HUGE); 389936Ssam } 399936Ssam x = frexp(arg,&exp); 409936Ssam while(x<0.5) { 419936Ssam x = x*2; 429936Ssam exp = exp-1; 439936Ssam } 449936Ssam if(x<sqrto2) { 459936Ssam x = 2*x; 469936Ssam exp = exp-1; 479936Ssam } 489936Ssam 499936Ssam z = (x-1)/(x+1); 509936Ssam zsq = z*z; 519936Ssam 529936Ssam temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0; 539936Ssam temp = temp/(((1.0*zsq + q2)*zsq + q1)*zsq + q0); 549936Ssam temp = temp*z + exp*log2; 559936Ssam return(temp); 569936Ssam } 579936Ssam 589936Ssam double log10(arg)599936Ssamlog10(arg) 609936Ssam double arg; 619936Ssam { 629936Ssam 639936Ssam return(log(arg)/ln10); 649936Ssam } 65