#StackBounty: #java #performance #reinventing-the-wheel #library #numerical-methods Arbitrary Precision nth Principal Root in Java – Ma…

Bounty: 50

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.
 * <p>
 * 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 * (31 - Integer.numberOfLeadingZeros(newPrecision)) >>> 1;

        // Make an initial guess:
        BigDecimal x = guessRoot(a, n);
        BigDecimal x0;

        // 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);
        }

        return x.round(c);
    }

    /**
     * Constructs an initial guess for the n-th principal root of
     * a given positive value.
     *
     * @param a The value whose n-th root is sought.
     * @param n The root to find.
     * @return An initial guess.
     */
    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);
    }

    /**
     * Returns the principal n-th root of the given positive value.
     *
     * @param decimal The value whose n-th root is sought.
     * @param n       The value of n needed.
     * @param context The MathContext to specify the precision and RoundingMode.
     * @return The principal n-th root of {@code decimal}.
     * @throws ArithmeticException If n is lesser than 2 or {@code decimal} is negative.
     */
    static BigDecimal principalRoot(final BigDecimal decimal,
                                    final int n,
                                    final MathContext context)
            throws ArithmeticException {
        if (n < 2)
            throw new ArithmeticException("'n' must at least be 2.");
        // Quick exits
        if (decimal.signum() == 0)
            return ZERO;
        if (decimal.compareTo(ONE) == 0)
            return ONE;
        if (decimal.signum() < 0)
            throw new ArithmeticException("Only positive values are supported.");
        return nthRoot(decimal, n, context);
    }
}

The method actually used is the last one: static BigDecimal principalRoot(...) in BigMath.java

Relevant methods of BigMath.java

... truncated ...

/**

 * Returns the square root of the given positive {@link BigDecimal} value.
 * The result has two extra bits of precision to ensure better accuracy.
 *
 * @param decimal The value whose square root is sought.
 * @param context The MathContext to specify the precision and RoundingMode.
 * @return The square root of {@code decimal}.
 * @throws ArithmeticException If {@code decimal} is negative.
 */
public static BigDecimal sqrt(final BigDecimal decimal,
                              final MathContext context)
        throws ArithmeticException {
    return Roots.principalRoot(decimal, 2, context);
}

/**
 * Returns the principal n-th root of the given positive value.
 *
 * @param decimal The value whose n-th root is sought.
 * @param n       The value of n needed.
 * @param context The MathContext to specify the precision and RoundingMode.
 * @return The principal n-th root of {@code decimal}.
 * @throws ArithmeticException If n is lesser than 2 or {@code decimal} is negative.
 */
public static BigDecimal principalRoot(final BigDecimal decimal,
                                       final int n,
                                       final MathContext context) {
    return Roots.principalRoot(decimal, n, context);
}

/**
 * A utility method that helps obtain a new {@link MathContext} from an existing
 * one. The new Context has the new precision specified but retains the {@link java.math.RoundingMode}.
 * <p>
 * Usually, it is used to increase the precision and hence "expand" the Context.
 *
 * @param c0           The initial {@link MathContext}.
 * @param newPrecision The required precision.
 * @return The new expanded Context.
 */
public static MathContext expandContext(MathContext c0, int newPrecision) {
    return new MathContext(
            newPrecision,
            c0.getRoundingMode()    // Retain rounding mode
    );
}
... truncated ...

Here’s sample usage with output

public static void main(String[] args) {
    int digits = 100;
    BigDecimal two = BigDecimal.valueOf(2);
    MathContext context = new MathContext(digits, RoundingMode.HALF_EVEN);

    BigDecimal root2 = BigMath.sqrt(two, context);
    System.out.println(root2.round(context));

    BigDecimal square = root2.pow(2, context);
    System.out.println(square);
}

Output

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573
2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

I welcome comments on all aspects of the code. I’d especially appreciate points on any cases which I’ve overlooked or ways to break my code.

My only request is please do not re-iterate the fact that I’m reinventing the wheel. I’m doing it so that I can learn.


Get this bounty!!!