home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / jpl_c / cbrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-05-18  |  1.6 KB  |  56 lines
  1. /* 1.0  04-27-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************/
  7.  
  8. #include "defs.h"
  9. #include "stdtyp.h"
  10. #include "errno.h"
  11. #include "mathcons.h"
  12.  
  13. /*
  14. Reference: See J. F. Hart, E. W. Cheney, C. L. Lawson, H. J. Maehly,
  15. C. K. Mesztenyi, J. R. Rice, H. C. Thacher, Jr., and C. Witzgall,
  16. COMPUTER APPROXIMATIONS, John Wiley & Sons, Inc., New York, NY, 1968.
  17. */
  18.  
  19. /************************************************************************/
  20.     double
  21. cbrt(x)            /* return cube root of x            */
  22.  
  23. /*----------------------------------------------------------------------*/
  24. double x;
  25. {
  26.     double fx, yn, frexp(), ldexp();
  27.     int n, ex, sgn;
  28.     
  29.     if (x IS 0.0)
  30.         return x;
  31.  
  32.     sgn = (x < 0.0 ? 1 : 0);
  33.     if (sgn)
  34.         x = -x;
  35.     fx = frexp(x, &ex);        /* frac x and exp x        */
  36.     yn = 0.59130052 + 0.41533799 * fx;
  37.         /* See Chapter 7, Table 640 of Reference. Has err=2.1 D */
  38.     yn = ldexp((yn + 1.5 * fx / (yn * yn + ldexp(fx, -1) / yn)), -1);
  39.         /* See Newton iteration formula 6.1.13 of Reference,    */
  40.     yn -= ldexp((yn - 1.5 * fx / (yn * yn + ldexp(fx, -1) / yn)), -1);
  41.         /* after two iterations, err = 19 D, without roundoff    */
  42.  
  43.     n = (ex % 3);
  44.     ex -= n;        /* ex now divisible by 3        */
  45.     if (n < 0)        /* must compensate if negative        */
  46.     {    n += 3;
  47.         ex -= 3;
  48.     }
  49.     if (n IS 1)
  50.         yn *= CBRTTWO;
  51.     else if (n IS 2)
  52.         yn *= CBRTFOUR;
  53.     yn = ldexp(yn, ex/3);
  54.     return (sgn ? -yn : yn);
  55. }
  56.