6
\$\begingroup\$

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.

\$\endgroup\$
2
+50
\$\begingroup\$
     * 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.

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ (+1) Points well worth consideration. I will definitely look into the secant method and improve the root guessing and comments. The reason I used the for loop and not the while (as you suggest) is that I had used it before to find several cases where the delta would not decrease any further and the program would run forever. The iteration limit ensures that my program halts. \$\endgroup\$ – Hungry Blue Dev Feb 24 '17 at 11:54
  • \$\begingroup\$ The secant method seems to have superlinear but subquadratic convergence (\$\phi\$). Does it really outperform Newton-Raphson in practical contexts for square roots? \$\endgroup\$ – Hungry Blue Dev Feb 24 '17 at 12:02
  • 1
    \$\begingroup\$ @Astrobleme, I can believe that the delta would not reach zero, but it doesn't need to reach zero to break the loop. It just needs to be less than 0.5 ULP at c0. \$\endgroup\$ – Peter Taylor Feb 24 '17 at 12:03
  • \$\begingroup\$ Oh. I was expecting delta to get smaller than 0.1 ULP at c0. (-‸ლ) \$\endgroup\$ – Hungry Blue Dev Feb 24 '17 at 12:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.