Arbitrary Precision nth Principal Root in Java - MathCore #1

This post is the first in the MathCore series.

The next post is here: Arbitrary precision π (Circular Constant) in Java - MathCore #2

Disclaimer

My project is too big to be reviewed in a single question. So, I'll post one class at a time.

The relevant methods needed from other classes are added as snippets.

Roots.java

The purpose of this class is to calculate the principal n-th root of a non-negative BigDecimal with the specified precision.

The precision used internally is actually a bit more so that the answer returned is fully accurate up to the specified precision and rounding can be done with confidence.

```
package mathcore.ops;

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;

import static java.math.BigDecimal.ONE;
import static java.math.BigDecimal.ZERO;
import static java.math.BigDecimal.valueOf;

/**
* A portion of BigMath refactored out to reduce overall complexity.
*
* This class handles the calculation of n-th principal roots.
*
* @author Subhomoy Haldar
* @version 1.0
*/
class Roots {
// Make this class un-instantiable
private Roots() {
}

/**
* Uses the n-th root algorithm to find principal root of a verified value.
*
* @param a The value whose n-th root is sought.
* @param n The root to find.
* @param c0 The initial (unexpanded) MathContext.
* @return The required principal root.
*/
private static BigDecimal nthRoot(final BigDecimal a,
final int n,
final MathContext c0) {
// Obtain constants that will be used in every iteration
final BigDecimal N = valueOf(n);
final int n_1 = n - 1;

// Increase precision by "n";
final int newPrecision = c0.getPrecision() + n;

final MathContext c = BigMath.expandContext(c0, newPrecision);

// The iteration limit (quadratic convergence)
final int limit = n n (3

* Uses the n-th root algorithm to find principal root of a verified value.

There's more than one algorithm to find an nth root. It would be more accurate to say "Uses the Newton-Raphson algorithm..."

I wonder whether it would be worth using the secant method instead, but I'm just throwing that out as an idea rather than a recommendation.

// Increase precision by "n";
        final int newPrecision = c0.getPrecision() + n;

This is a classic example of a bad comment. Anyone can see that the code increases precision by n, but what isn't obvious is why by n? I would think it makes sense to increase the precision by the logarithm of n to some base.

// The iteration limit (quadratic convergence)
        final int limit = n * n * (31 - Integer.numberOfLeadingZeros(newPrecision)) >>> 1;

This comment is better, but a reference or sketch proof for why this is the appropriate value would be good.

Also, I'm not convinced that this is the best stopping condition in general. See below.

// Iterate
        for (int i = 0; i < limit; i++) {
            x0 = x;
            BigDecimal delta = a.divide(x0.pow(n_1), c)
                    .subtract(x0, c)
                    .divide(N, c);
            x = x0.add(delta, c);
        }

A comment along the lines of

// Newton-Raphson to find zero of f(x) = x^n - a
        // x' = x - f(x) / f'(x) = x - (x^n - a) / (nx^{n-1})

would be nice.

I'm not convinced of the benefit of renaming x to x0 for a scope in which x doesn't change value.

Since you have delta, you could break out of the loop if you've converged before limit iterations. I would expect this to give a significant speedup.

return x.round(c);

Is that a bug? I think it should be x.round(c0).

private static BigDecimal guessRoot(BigDecimal a, int n) {
        // 1. Obtain first (1/n)th of total bits of magnitude
        BigInteger magnitude = a.unscaledValue();
        final int length = magnitude.bitLength() * (n - 1) / n;
        magnitude = magnitude.shiftRight(length);

        // 2. Obtain the correct scale
        final int newScale = a.scale() / n;

        // 3. Construct the initial guess
        return new BigDecimal(magnitude, newScale);
    }

That's one way of doing it. Have you considered the alternative of

double lowPrecisionA = a.doubleValue();
        double lowPrecisionRoot = Math.exp(Math.log(a) / n);
        return BigDecimal.valueOf(lowPrecisionRoot)

? It should get you about fourteen or fifteen more decimal digits of accurate initial value, saving three or four expensive iterations of Newton-Raphson.

StackExchange Code Review Q#151676, answer score: 2