Return to log.c CVS log

Up to [Research Unix] / researchv9 / libc / math

Annotation of researchv9/libc/math/log.c, revision 1.1.1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.