Subversion Repositories Kolibri OS

Rev

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

  1. /*                                                      cbrtl.c
  2.  *
  3.  *      Cube root, long double precision
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * long double x, y, cbrtl();
  10.  *
  11.  * y = cbrtl( 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.  *    IEEE     .125,8        80000      7.0e-20     2.2e-20
  33.  *    IEEE    exp(+-707)    100000      7.0e-20     2.4e-20
  34.  *
  35.  */
  36.  
  37. /*
  38. Cephes Math Library Release 2.2: January, 1991
  39. Copyright 1984, 1991 by Stephen L. Moshier
  40. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  41. */
  42.  
  43. /*
  44.   Modified for mingwex.a
  45.   2002-07-01  Danny Smith  <dannysmith@users.sourceforge.net>
  46.  */
  47. #ifdef __MINGW32__
  48. #include "cephes_mconf.h"
  49. #else
  50. #include "mconf.h"
  51. #endif
  52.  
  53. static const long double CBRT2  = 1.2599210498948731647672L;
  54. static const long double CBRT4  = 1.5874010519681994747517L;
  55. static const long double CBRT2I = 0.79370052598409973737585L;
  56. static const long double CBRT4I = 0.62996052494743658238361L;
  57.  
  58. #ifndef __MINGW32__
  59.  
  60. #ifdef ANSIPROT
  61. extern long double frexpl ( long double, int * );
  62. extern long double ldexpl ( long double, int );
  63. extern int isnanl ( long double );
  64. #else
  65. long double frexpl(), ldexpl();
  66. extern int isnanl();
  67. #endif
  68.  
  69. #ifdef INFINITIES
  70. extern long double INFINITYL;
  71. #endif
  72.  
  73. #endif /* __MINGW32__ */
  74.  
  75. long double cbrtl(x)
  76. long double x;
  77. {
  78. int e, rem, sign;
  79. long double z;
  80.  
  81. #ifdef __MINGW32__
  82. if (!isfinite (x) || x == 0.0L)
  83.         return(x);
  84. #else    
  85.  
  86. #ifdef NANS
  87. if(isnanl(x))
  88.         return(x);
  89. #endif
  90. #ifdef INFINITIES
  91. if( x == INFINITYL)
  92.         return(x);
  93. if( x == -INFINITYL)
  94.         return(x);
  95. #endif
  96. if( x == 0 )
  97.         return( x );
  98.  
  99. #endif /* __MINGW32__ */
  100.  
  101. if( x > 0 )
  102.         sign = 1;
  103. else
  104.         {
  105.         sign = -1;
  106.         x = -x;
  107.         }
  108.  
  109. z = x;
  110. /* extract power of 2, leaving
  111.  * mantissa between 0.5 and 1
  112.  */
  113. x = frexpl( x, &e );
  114.  
  115. /* Approximate cube root of number between .5 and 1,
  116.  * peak relative error = 1.2e-6
  117.  */
  118. x = (((( 1.3584464340920900529734e-1L * x
  119.        - 6.3986917220457538402318e-1L) * x
  120.        + 1.2875551670318751538055e0L) * x
  121.        - 1.4897083391357284957891e0L) * x
  122.        + 1.3304961236013647092521e0L) * x
  123.        + 3.7568280825958912391243e-1L;
  124.  
  125. /* exponent divided by 3 */
  126. if( e >= 0 )
  127.         {
  128.         rem = e;
  129.         e /= 3;
  130.         rem -= 3*e;
  131.         if( rem == 1 )
  132.                 x *= CBRT2;
  133.         else if( rem == 2 )
  134.                 x *= CBRT4;
  135.         }
  136. else
  137.         { /* argument less than 1 */
  138.         e = -e;
  139.         rem = e;
  140.         e /= 3;
  141.         rem -= 3*e;
  142.         if( rem == 1 )
  143.                 x *= CBRT2I;
  144.         else if( rem == 2 )
  145.                 x *= CBRT4I;
  146.         e = -e;
  147.         }
  148.  
  149. /* multiply by power of 2 */
  150. x = ldexpl( x, e );
  151.  
  152. /* Newton iteration */
  153.  
  154. x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
  155. x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
  156.  
  157. if( sign < 0 )
  158.         x = -x;
  159. return(x);
  160. }
  161.