ph-math/src/main/java/com/helger/numbercruncher/mathutils/BigFunctions.java
public static BigDecimal intRoot (@Nonnull final BigDecimal px, final long index, final int scale)
  {
    BigDecimal x = px;

    // Check that x >= 0.
    if (x.signum () < 0)
    {
      throw new IllegalArgumentException ("x < 0: " + x);
    }

    final int sp1 = scale + 1;
    final BigDecimal n = x;
    final BigDecimal i = BigDecimal.valueOf (index);
    final BigDecimal im1 = BigDecimal.valueOf (index - 1);
    final BigDecimal tolerance = BigDecimal.valueOf (5).movePointLeft (sp1);
    BigDecimal xPrev;

    // The initial approximation is x/index.
    x = x.divide (i, scale, RoundingMode.HALF_EVEN);

    // Loop until the approximations converge
    // (two successive approximations are equal after rounding).
    do
    {
      // x^(index-1)
      final BigDecimal xToIm1 = intPower (x, index - 1, sp1);

      // x^index
      final BigDecimal xToI = x.multiply (xToIm1).setScale (sp1, RoundingMode.HALF_EVEN);

      // n + (index-1)*(x^index)
      final BigDecimal numerator = n.add (im1.multiply (xToI)).setScale (sp1, RoundingMode.HALF_EVEN);

      // (index*(x^(index-1))
      final BigDecimal denominator = i.multiply (xToIm1).setScale (sp1, RoundingMode.HALF_EVEN);

      // x = (n + (index-1)*(x^index)) / (index*(x^(index-1)))
      xPrev = x;
      x = numerator.divide (denominator, sp1, RoundingMode.DOWN);

      Thread.yield ();
    } while (x.subtract (xPrev).abs ().compareTo (tolerance) > 0);

    return x;
  }
public static BigDecimal intRoot (@Nonnull final BigDecimal px, final long index, final int scale)
  {
    BigDecimal x = px;

    // Check that x >= 0.
    if (x.signum () < 0)
    {
      throw new IllegalArgumentException ("x < 0: " + x);
    }

    final int sp1 = scale + 1;
    final BigDecimal n = x;
    final BigDecimal i = BigDecimal.valueOf (index);
    final BigDecimal im1 = BigDecimal.valueOf (index - 1);
    final BigDecimal tolerance = BigDecimal.valueOf (5).movePointLeft (sp1);
    BigDecimal xPrev;

    // The initial approximation is x/index.
    x = x.divide (i, scale, RoundingMode.HALF_EVEN);

    // Loop until the approximations converge
    // (two successive approximations are equal after rounding).
    do
    {
      // x^(index-1)
      final BigDecimal xToIm1 = intPower (x, index - 1, sp1);

      // x^index
      final BigDecimal xToI = x.multiply (xToIm1).setScale (sp1, RoundingMode.HALF_EVEN);

      // n + (index-1)*(x^index)
      final BigDecimal numerator = n.add (im1.multiply (xToI)).setScale (sp1, RoundingMode.HALF_EVEN);

      // (index*(x^(index-1))
      final BigDecimal denominator = i.multiply (xToIm1).setScale (sp1, RoundingMode.HALF_EVEN);

      // x = (n + (index-1)*(x^index)) / (index*(x^(index-1)))
      xPrev = x;
      x = numerator.divide (denominator, sp1, RoundingMode.DOWN);

      Thread.yield ();
    } while (x.subtract (xPrev).abs ().compareTo (tolerance) > 0);

    return x;
  }
Compute the integral root of x to a given scale, x &ge; 0. Use Newton's
algorithm.

@param px
the value of x
@param index
the integral root value
@param scale
the desired scale of the result
@return the result value