/*
   ###########################################################################
   This function computes the arbitrary-precision cube root of a given
   numerical string X, rounded to any specified number of decimals up
   to a maximum of 100 decimals.

   CAVEAT:
   There may be a few +/- fuzzy digits at the end of extremely long values.

   ARGUMENTS:
   X = Numerical string for which we seek the cube root.

   Decimals = Number of decimals at which to round off the result
              if exact resolution is not possible.

   This function is based on the following cube root iteration algorithm.

   -------------------------
   CUBE ROOT of X ITERATION:

   Let:
   a = Current approximation to cube root of (X).
       Initial approximation = Normal native machine-precision value.

   b = Next generation approximation computed from current (a).

   p = Precision limit.
       10^(-(decimals + 1))
       If we want 16 decimals of precision,
       then we set:
       p = 10^(-(16 + 1))
         = 10^(-17)
         = 10E+17
         = 0.00000000000000001

   At the end of each iteration, the resulting value of (b) becomes the
   new value of (a) for the next cycle and the process is repeated until
   the desired level of accuracy is reached.

   The iteration process continues while (abs(b-a) >= p)

   When the difference abs(b-a) < p, then we are
   done and the current value of (b) is returned.

   -----------------------------------------------------------
   CUBE ROOT of X ITERATION ALGORITHM FOR ARBITRARY-PRECISION

   The initial expression:
   Approx = X^(1/3)
   refers to the standard 16-digit double-precision value used
   as the first approximation to the arbitrary-precision value
   of the cube root of X.
   -----------------------------------------------------------
   ARGUMENTS (X, Decimals)

  _START_

   Approx = X^(1/3)
   Diff = 1
   PrecisionLimit = 10^(-(Decimals + 1))

   LOOP_WHILE (Diff >= PrecisionLimit)
  {
   CubeRoot = (Approx + X/(Approx^2)) / 2
   Diff = abs(CubeRoot - Approx)
   Approx = CubeRoot
  }
   END_LOOP

   return CubeRoot

  _END_


   ERRORS:
   FALSE is returned if either argument is non-numeric.

   NO DEPENDENCIES
   ###########################################################################
*/

   function BC_Cube_Root ($XStr, $Decimals=16)
{
   $X = trim($XStr);

   $Decimals = trim($Decimals);

// ----------------------------------------
// Error if either argument is non-numeric.

   if (!Is_Numeric($X) or !Is_Numeric($Decimals))  {return FALSE;}

// -------------------------------------------
// Error if Decimals argument is < 1 or > 100.

   if ($Decimals < 1 or $Decimals > 100)
      {
       return FALSE;
      }

   $Decimals = floor(abs($Decimals));

// ------------------------------------------------------------
// Generate 1st approximation = Native machine-precision value.

   $Approx = SPrintF("%1.16f", pow($X, 1/3));
   $D = 100;

/* ----------------------------------------------------------
   Loop to iterate cube root value (limit = 100 iterations).
   This should never get anywhere near 100.  This limit helps
   to prevent an infinite loop if something goes wrong.
*/
   for ($i=1;  $i <= 100;   $i++)
  {
// --------------------------------------------------------------
// Compute next approximation (b) from current approximation (a).

   $CubeRoot = bcDiv(bcAdd(bcMul('2', $Approx, $D),
               bcDiv($X, bcPow($Approx, '2', $D),$D),$D), '3', $D);

/* --------------------------------------------------------------
   If (Approx == CubeRoot), to the set precision limit, then done.
   Otherwise, update Approx and CubeRoot and execute another cycle.
*/
   if (bcComp($Approx,$CubeRoot,$D) == 0){break;}
   else
      {$Approx = $CubeRoot;}
  }

/* ------------------------------------------
   Done.  Return the cube root value rounded
   to the specified number of decimals.
*/
   return RTrim(RTrim(bcAdd($CubeRoot, '0.' . Str_Repeat(0, $Decimals)
   . '5', $Decimals), '0'), '.');

} // End of  BC_Cube_Root (...)