Java - MathCore #1中任意精度的第n主根

Question

这篇文章是MathCore系列中的第一篇.

下一篇文章是：( MathCore #2中的任意精度π(圆形常数)

免责声明

我的项目太大了，不能用一个问题来复习。所以，我会一次发一个班。

从其他类中需要的相关方法作为代码片段添加。

Roots.java

此类的目的是计算具有指定精度的非负BigDecimal的主n根。

内部使用的精度实际上要高一点，因此返回的答案完全准确，直到指定的精度，并且可以有信心地进行舍入。

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

实际使用的方法是最后一个方法：static BigDecimal principalRoot(...) in BigMath.java

相关的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 ...

下面是输出的示例用法

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

输出

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573
2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

我欢迎就守则的所有方面发表意见。我特别想知道我忽略的任何情况或破坏我的代码的方法。

我唯一的要求是，请不要重复我正在重新发明车轮的事实。我这么做是为了让我学会。

Answer 1

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

找到第n根的算法不止一种。更准确的说法是“使用牛顿-拉夫森算法.”

我不知道它是否值得使用割线方法，但我只是把它作为一个想法而不是一个建议。

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

这是一个典型的坏评论的例子。任何人都可以看到代码增加了n的精度，但不明显的是为什么增加了n？我认为用n的对数将精度提高到某个基数是有意义的。

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

这个评论是更好的，但参考或草图证明为什么这是适当的价值将是很好的。

而且，我不相信这是最好的停止条件在一般情况下。见下文。

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

与…有关的评论

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

会很好的。

我不相信将x重命名为x0的好处，因为在这个范围中x不改变值。

由于您有delta，所以如果您在limit迭代之前进行了收敛，就可以脱离循环。我希望这能大大加快速度。

        return x.round(c);

那是窃听器吗？我想应该是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);
    }

这是一种方法。你有没有考虑过

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

？它会让你多得到十四到十五位精确初值的十进制数字，省去了牛顿-拉夫森的三到四次昂贵的迭代。