|
|
1.1 ! root 1: /* ! 2: * log - calculate the natural logarithm of arg. ! 3: * ! 4: * The coefficients are #2705 from Hart & Cheney. (19.38D) ! 5: */ ! 6: ! 7: #include <errno.h> ! 8: #include <math.h> ! 9: ! 10: int errno; ! 11: double frexp(); ! 12: static double log2 = 0.693147180559945309e0; ! 13: static double ln10 = 2.302585092994045684; ! 14: static double sqrto2 = 0.707106781186547524e0; ! 15: static double p0 = -.240139179559210510e2; ! 16: static double p1 = 0.309572928215376501e2; ! 17: static double p2 = -.963769093368686593e1; ! 18: static double p3 = 0.421087371217979714e0; ! 19: static double q0 = -.120069589779605255e2; ! 20: static double q1 = 0.194809660700889731e2; ! 21: static double q2 = -.891110902798312337e1; ! 22: ! 23: double ! 24: log(arg) ! 25: double arg; ! 26: { ! 27: double x, z, zsq, temp; ! 28: int exp; ! 29: ! 30: if(arg <= 0.) { ! 31: errno = EDOM; ! 32: return(-HUGE); ! 33: } ! 34: x = frexp(arg,&exp); ! 35: while(x<0.5) { ! 36: x = x*2; ! 37: exp = exp-1; ! 38: } ! 39: if(x<sqrto2) { ! 40: x = 2*x; ! 41: exp = exp-1; ! 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/(((1.0*zsq + q2)*zsq + q1)*zsq + q0); ! 49: temp = temp*z + exp*log2; ! 50: return(temp); ! 51: } ! 52: ! 53: double ! 54: log10(arg) ! 55: double arg; ! 56: { ! 57: ! 58: return(log(arg)/ln10); ! 59: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.