sys/src/libc/port/log.c

   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 <u.h>
  11 #include <libc.h>
  12
  13 #define log2    0.693147180559945309e0
  14 #define ln10o1  .4342944819032518276511
  15 #define sqrto2  0.707106781186547524e0
  16 #define p0      -.240139179559210510e2
  17 #define p1      0.309572928215376501e2
  18 #define p2      -.963769093377840513e1
  19 #define p3      0.421087371217979714e0
  20 #define q0      -.120069589779605255e2
  21 #define q1      0.194809660700889731e2
  22 #define q2      -.891110902798312337e1
  23
  24 double
  25 log(double arg)
  26 {
  27         double x, z, zsq, temp;
  28         int exp;
  29
  30         if(arg <= 0)
  31                 return NaN();
  32         x = frexp(arg, &exp);
  33         while(x < 0.5) {
  34                 x *= 2;
  35                 exp--;
  36         }
  37         if(x < sqrto2) {
  38                 x *= 2;
  39                 exp--;
  40         }
  41
  42         z = (x-1) / (x+1);
  43         zsq = z*z;
  44
  45         temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
  46         temp = temp/(((zsq + q2)*zsq + q1)*zsq + q0);
  47         temp = temp*z + exp*log2;
  48         return temp;
  49 }
  50
  51 double
  52 log10(double arg)
  53 {
  54
  55         if(arg <= 0)
  56                 return NaN();
  57         return log(arg) * ln10o1;
  58 }