Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  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        -10,10     200000      1.8e-17     6.2e-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.2:  January, 1991
  40. Copyright 1984, 1991 by Stephen L. Moshier
  41. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  42. */
  43.  
  44. /*
  45.   Modified for mingwex.a
  46.   2002-07-01  Danny Smith  <dannysmith@users.sourceforge.net>
  47.  */
  48. #ifdef __MINGW32__
  49. #include <math.h>
  50. #include "cephes_mconf.h"
  51. #else
  52. #include "mconf.h"
  53. #endif
  54.  
  55.  
  56. static const double CBRT2  = 1.2599210498948731647672;
  57. static const double CBRT4  = 1.5874010519681994747517;
  58. static const double CBRT2I = 0.79370052598409973737585;
  59. static const double CBRT4I = 0.62996052494743658238361;
  60.  
  61. #ifndef __MINGW32__
  62. #ifdef ANSIPROT
  63. extern double frexp ( double, int * );
  64. extern double ldexp ( double, int );
  65. extern int isnan ( double );
  66. extern int isfinite ( double );
  67. #else
  68. double frexp(), ldexp();
  69. int isnan(), isfinite();
  70. #endif
  71. #endif
  72.  
  73. double cbrt(x)
  74. double x;
  75. {
  76. int e, rem, sign;
  77. double z;
  78.  
  79. #ifdef __MINGW32__
  80. if (!isfinite (x) || x == 0 )
  81.   return x;
  82. #else
  83.  
  84. #ifdef NANS
  85. if( isnan(x) )
  86.   return x;
  87. #endif
  88. #ifdef INFINITIES
  89. if( !isfinite(x) )
  90.   return x;
  91. #endif
  92. if( x == 0 )
  93.         return( x );
  94.  
  95. #endif /* __MINGW32__ */
  96.  
  97. if( x > 0 )
  98.         sign = 1;
  99. else
  100.         {
  101.         sign = -1;
  102.         x = -x;
  103.         }
  104.  
  105. z = x;
  106. /* extract power of 2, leaving
  107.  * mantissa between 0.5 and 1
  108.  */
  109. x = frexp( x, &e );
  110.  
  111. /* Approximate cube root of number between .5 and 1,
  112.  * peak relative error = 9.2e-6
  113.  */
  114. x = (((-1.3466110473359520655053e-1  * x
  115.       + 5.4664601366395524503440e-1) * x
  116.       - 9.5438224771509446525043e-1) * x
  117.       + 1.1399983354717293273738e0 ) * x
  118.       + 4.0238979564544752126924e-1;
  119.  
  120. /* exponent divided by 3 */
  121. if( e >= 0 )
  122.         {
  123.         rem = e;
  124.         e /= 3;
  125.         rem -= 3*e;
  126.         if( rem == 1 )
  127.                 x *= CBRT2;
  128.         else if( rem == 2 )
  129.                 x *= CBRT4;
  130.         }
  131.  
  132.  
  133. /* argument less than 1 */
  134.  
  135. else
  136.         {
  137.         e = -e;
  138.         rem = e;
  139.         e /= 3;
  140.         rem -= 3*e;
  141.         if( rem == 1 )
  142.                 x *= CBRT2I;
  143.         else if( rem == 2 )
  144.                 x *= CBRT4I;
  145.         e = -e;
  146.         }
  147.  
  148. /* multiply by power of 2 */
  149. x = ldexp( x, e );
  150.  
  151. /* Newton iteration */
  152. x -= ( x - (z/(x*x)) )*0.33333333333333333333;
  153. #ifdef DEC
  154. x -= ( x - (z/(x*x)) )/3.0;
  155. #else
  156. x -= ( x - (z/(x*x)) )*0.33333333333333333333;
  157. #endif
  158.  
  159. if( sign < 0 )
  160.         x = -x;
  161. return(x);
  162. }
  163.