home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cmath / cbrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.5 KB  |  148 lines
  1. /*                            cbrt.c
  2.  *
  3.  *    Cube root
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, cbrt();
  10.  *
  11.  * y = cbrt( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the cube root of the argument, which may be negative.
  18.  *
  19.  * Range reduction involves determining the power of 2 of
  20.  * the argument.  A polynomial of degree 2 applied to the
  21.  * mantissa, and multiplication by the cube root of 1, 2, or 4
  22.  * approximates the root to within about 0.1%.  Then Newton's
  23.  * iteration is used three times to converge to an accurate
  24.  * result.
  25.  *
  26.  *
  27.  *
  28.  * ACCURACY:
  29.  *
  30.  *                      Relative error:
  31.  * arithmetic   domain     # trials      peak         rms
  32.  *    DEC        0,8        50000       1.8e-17     5.6e-18
  33.  *    IEEE       0,1e308    30000       1.5e-16     5.0e-17
  34.  *
  35.  */
  36. /*                            cbrt.c  */
  37.  
  38. /*
  39. Cephes Math Library Release 2.1:  December, 1988
  40. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  41. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  42. */
  43.  
  44.  
  45. #include "mconf.h"
  46.  
  47. #ifdef UNK
  48. static double CBRT2 = 1.25992104989487316477;
  49. static double CBRT4 = 1.58740105196819947475;
  50. #endif
  51.  
  52. #ifdef DEC
  53. static short CBT2[] = {040241,042427,0146153,0112127};
  54. static short CBT4[] = {040313,027765,024753,070744};
  55. #define CBRT2 *(double *)CBT2
  56. #define CBRT4 *(double *)CBT4
  57. #endif
  58.  
  59. #ifdef IBMPC
  60. static short CBT2[] = {0x728b,0xf98d,0x28a2,0x3ff4};
  61. static short CBT4[] = {0x6e3d,0xa53d,0x65fe,0x3ff9};
  62. #define CBRT2 *(double *)CBT2
  63. #define CBRT4 *(double *)CBT4
  64. #endif
  65.  
  66. #ifdef MIEEE
  67. static short CBT2[] = {
  68. 0x3ff4,0x28a2,0xf98d,0x728b
  69. };
  70. static short CBT4[] = {
  71. 0x3ff9,0x65fe,0xa53d,0x6e3d,
  72. };
  73. #define CBRT2 *(double *)CBT2
  74. #define CBRT4 *(double *)CBT4
  75. #endif
  76.  
  77. double cbrt(x)
  78. double x;
  79. {
  80. int e, rem, sign;
  81. double z;
  82. double frexp(), ldexp();
  83.  
  84. if( x == 0 )
  85.     return( 0.0 );
  86. if( x > 0 )
  87.     sign = 1;
  88. else
  89.     {
  90.     sign = -1;
  91.     x = -x;
  92.     }
  93.  
  94. z = x;
  95. /* extract power of 2, leaving
  96.  * mantissa between 0.5 and 1
  97.  */
  98. x = frexp( x, &e );
  99.  
  100. /* Approximate cube root of number between .5 and 1,
  101.  * peak relative error = 6.36e-4
  102.  */
  103. x = (-1.9150215751434832257e-1  * x
  104.     + 6.9757045195246484393e-1) * x
  105.     + 4.9329566506409572486e-1;
  106.  
  107.  
  108. /* exponent divided by 3 */
  109. if( e >= 0 )
  110.     {
  111.     rem = e;
  112.     e /= 3;
  113.     rem -= 3*e;
  114.     if( rem == 1 )
  115.         x *= CBRT2;
  116.     else if( rem == 2 )
  117.         x *= CBRT4;
  118.     }
  119.  
  120.  
  121. /* argument less than 1 */
  122.  
  123. else
  124.     {
  125.     e = -e;
  126.     rem = e;
  127.     e /= 3;
  128.     rem -= 3*e;
  129.     if( rem == 1 )
  130.         x /= CBRT2;
  131.     else if( rem == 2 )
  132.         x /= CBRT4;
  133.     e = -e;
  134.     }
  135.  
  136. /* multiply by power of 2 */
  137. x = ldexp( x, e );
  138.  
  139. /* Newton iteration */
  140. x -= ( x - (z/(x*x)) )/3.0;
  141. x -= ( x - (z/(x*x)) )/3.0;
  142. x -= ( x - (z/(x*x)) )/3.0;
  143.  
  144. if( sign < 0 )
  145.     x = -x;
  146. return(x);
  147. }
  148.