Return to log.c CVS log

Up to [CSRG BSD Unix] / 43BSD / contrib / icon / rt

Annotation of 43BSD/contrib/icon/rt/log.c, revision 1.1.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:    }

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