home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / jpl_c / log.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-05-18  |  1.9 KB  |  70 lines
  1. /* 1.0  04-27-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************
  7.  *    Programmmed using the algorithms given in:
  8.  *
  9.  *    Coty, William J., Jr., and Waite, William, SOFTWARE MANUAL FOR
  10.  *    THE ELEMENTARY FUNCTIONS, Prentice-Hall Series in Computational
  11.  *    Mathematics, Prentice-Hall, Inc., Inglewood Cliffs, NJ, 1980,
  12.  *    pp. 35-59.
  13.  *
  14.  *----------------------------------------------------------------------*/
  15.  
  16. #include "defs.h"
  17. #include "stdtyp.h"
  18. #include "errno.h"
  19. #include "mathtyp.h"
  20. #include "mathcons.h"
  21.  
  22. /*----------------------------------------------------------------------*/
  23.  
  24. #define A0     -0.64124943423745581147e+2
  25. #define A1     +0.16383943563021534222e+2
  26. #define A2     -0.78956112887491257267e+0
  27. #define A(w)     ((A2*w A1)*w A0)
  28.  
  29. #define B0     -0.76949932108494879777e+3
  30. #define B1     +0.31203222091924532844e+3
  31. #define B2     -0.35667977739034646171e+2
  32. #define B3    +1.0
  33. #define B(w)    (((w B2)*w B1)*w B0)
  34.  
  35. #define C1     0.693359375            /* = 355 / 512        */
  36. #define C2     (-2.121944400546905827679e-4)    /* C1 + C2 = log(2)    */
  37.  
  38. /*\p*********************************************************************/
  39.     double
  40. log(x)            /* return natural log (base e) of x        */
  41.  
  42. /*----------------------------------------------------------------------*/
  43. double x;
  44. {
  45.     double r, f, z, w, znum, zden, xn;
  46.     int n;
  47.     
  48.     if (x <= 0.0)
  49.     {    errno = EDOM;
  50.         return -INFINITY;
  51.     }
  52.     f = frexp(x, &n);
  53.     if (f > ROOTHALF)
  54.     {    znum = (f - 0.5) - 0.5;
  55.         zden = ldexp(f, -1) + 0.5;
  56.     }
  57.     else
  58.     {    --n;
  59.         znum = f - 0.5;
  60.         zden = ldexp(znum, -1) + 0.5;
  61.     }
  62.     z = znum / zden;
  63.     w = z * z;
  64.     r = z + z * (w *
  65.          A(w)
  66.          /B(w));
  67.     xn = n;
  68.     return (xn * C2 + r) + xn * C1;
  69. }
  70.