home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / jpl_c / exp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-05-18  |  1.9 KB  |  71 lines
  1. /* 1.0  01-09-85 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1985        *
  6.  ************************************************************************
  7.  *    Programmmed using the algorithm 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. 60-83.
  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 P0     0.249999999999999993e+0
  25. #define P1     0.694360001511792852e-2
  26. #define P2     0.165203300268279130e-4
  27. #define P(z)    ((P2*z+P1)*z+P0)
  28.  
  29. #define Q0     0.5
  30. #define Q1     0.555538666969001188e-1
  31. #define Q2     0.495862884905441294e-3
  32. #define Q(z)     ((Q2*z+Q1)*z+Q0)
  33.  
  34. #define C1    0.693359375            /* = 355 / 512        */
  35. #define C2    2.121944400546905827697e-4    /* C1 - C2 = LOG2    */
  36.  
  37. /*\p*********************************************************************/
  38.     double
  39. exp(x)                /* exponential function of x, e^x    */
  40.  
  41. /*----------------------------------------------------------------------*/
  42. double x;
  43. {
  44.     double x2, xn, g, r, z;
  45.     int n;
  46.     
  47.     if (x > LOGINFINITY)
  48.     {    errno = ERANGE;
  49.         return INFINITY;
  50.     }
  51.     if (x < LOGLEAST)
  52.     {    errno = ERANGE;
  53.         return 0.0;
  54.     }
  55.     if (ABS(x) < WASHOUT)
  56.         return 1.0;
  57.  
  58.     x2 = modf(x * LOG2e, &xn);
  59.     if (x2 >= 0.5)
  60.         ++xn;
  61.     else if (x2 <= -0.5)
  62.         --xn;
  63.     n = xn;
  64.     x2 = modf(x, &x);        /* in case no guard digit (p. 68) */
  65.     g = ((x - xn * C1) + x2) + xn * C2;
  66.     z = g * g;
  67.     r = P(z) * g;
  68.     r = 0.5 + r / (Q(z) - r);
  69.     return ldexp(r, n + 1);
  70. }
  71.