lib/ack/libp/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
   8 /* $Header$ */
   9
  10 #define __NO_DEFS
  11 #include <math.h>
  12 #include <pc_err.h>
  13
  14 #if __STDC__
  15 #include <pc_math.h>
  16 #include <float.h>
  17 #endif
  18 #undef HUGE
  19 #define HUGE    1e1000
  20
  21 double
  22 _log(x)
  23         double  x;
  24 {
  25         /*      Algorithm and coefficients from:
  26                         "Software manual for the elementary functions"
  27                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  28         */
  29         static double a[] = {
  30                 -0.64124943423745581147e2,
  31                  0.16383943563021534222e2,
  32                 -0.78956112887491257267e0
  33         };
  34         static double b[] = {
  35                 -0.76949932108494879777e3,
  36                  0.31203222091924532844e3,
  37                 -0.35667977739034646171e2,
  38                  1.0
  39         };
  40
  41         extern double   _fef();
  42         double  znum, zden, z, w;
  43         int     exponent;
  44
  45         if (x <= 0) {
  46                 _trp(ELOG);
  47                 return -HUGE;
  48         }
  49
  50         x = _fef(x, &exponent);
  51         if (x > M_1_SQRT2) {
  52                 znum = (x - 0.5) - 0.5;
  53                 zden = x * 0.5 + 0.5;
  54         }
  55         else {
  56                 znum = x - 0.5;
  57                 zden = znum * 0.5 + 0.5;
  58                 exponent--;
  59         }
  60         z = znum/zden; w = z * z;
  61         x = z + z * w * (POLYNOM2(w,a)/POLYNOM3(w,b));
  62         z = exponent;
  63         x += z * (-2.121944400546905827679e-4);
  64         return x + z * 0.693359375;
  65 }