lib/math/log.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14 double
  15 log(double x)
  16 {
  17         /*      Algorithm and coefficients from:
  18                         "Software manual for the elementary functions"
  19                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  20         */
  21         static double a[] = {
  22                 -0.64124943423745581147e2,
  23                  0.16383943563021534222e2,
  24                 -0.78956112887491257267e0
  25         };
  26         static double b[] = {
  27                 -0.76949932108494879777e3,
  28                  0.31203222091924532844e3,
  29                 -0.35667977739034646171e2,
  30                  1.0
  31         };
  32
  33         double  znum, zden, z, w;
  34         int     exponent;
  35
  36         if (__IsNan(x)) {
  37                 errno = EDOM;
  38                 return x;
  39         }
  40         if (x < 0) {
  41                 errno = EDOM;
  42                 return -HUGE_VAL;
  43         }
  44         else if (x == 0) {
  45                 errno = ERANGE;
  46                 return -HUGE_VAL;
  47         }
  48
  49         if (x <= DBL_MAX) {
  50         }
  51         else return x;  /* for infinity and Nan */
  52         x = frexp(x, &exponent);
  53         if (x > M_1_SQRT2) {
  54                 znum = (x - 0.5) - 0.5;
  55                 zden = x * 0.5 + 0.5;
  56         }
  57         else {
  58                 znum = x - 0.5;
  59                 zden = znum * 0.5 + 0.5;
  60                 exponent--;
  61         }
  62         z = znum/zden; w = z * z;
  63         x = z + z * w * (POLYNOM2(w,a)/POLYNOM3(w,b));
  64         z = exponent;
  65         x += z * (-2.121944400546905827679e-4);
  66         return x + z * 0.693359375;
  67 }