Annotation of 43BSD/contrib/icon/rt/log.c, revision 1.1

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:    }

unix.superglobalmegacorp.com

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