Newton's method with specified digits of precision

Question

I'm trying to write a function in Java that calculates the n-th root of a number. I'm using Newton's method for this. However, the user should be able to specify how many digits of precision they want. This is the part with which I'm having trouble, as my answer is often not entirely correct. The relevant code is here: http://pastebin.com/d3rdpLW8. How could I fix this code so that it always gives the answer to at least p digits of precision? (without doing more work than is necessary)

import java.util.Random;

public final class Compute {

    private Compute() {
    }

    public static void main(String[] args) {
        Random rand = new Random(1230);
        for (int i = 0; i < 500000; i++) {
            double k = rand.nextDouble()/100;
            int n = (int)(rand.nextDouble() * 20) + 1;
            int p = (int)(rand.nextDouble() * 10) + 1;
            double math = n == 0 ? 1d : Math.pow(k, 1d / n);
            double compute = Compute.root(n, k, p);
            if(!String.format("%."+p+"f", math).equals(String.format("%."+p+"f", compute))) {
                System.out.println(String.format("%."+p+"f", math));
                System.out.println(String.format("%."+p+"f", compute));
                System.out.println(math + " " + compute + " " + p);
            }
        }
    }

    /**
     * Returns the n-th root of a positive double k, accurate to p decimal
     * digits.
     * 
     * @param n
     *            the degree of the root.
     * @param k
     *            the number to be rooted.
     * @param p
     *            the decimal digit precision.
     * @return the n-th root of k
     */
    public static double root(int n, double k, int p) {     
        double epsilon = pow(0.1, p+2);
        double approx = estimate_root(n, k);
        double approx_prev;

        do {
            approx_prev = approx;
            // f(x) / f'(x) = (x^n - k) / (n * x^(n-1)) = (x - k/x^(n-1)) / n
            approx -= (approx - k / pow(approx, n-1)) / n;
        } while (abs(approx - approx_prev) > epsilon);
        return approx;
    }

    private static double pow(double x, int y) {
        if (y == 0)
            return 1d;
        if (y == 1)
            return x;
        double k = pow(x * x, y >> 1);
        return (y & 1) == 0 ? k : k * x;
    }

    private static double abs(double x) {
        return Double.longBitsToDouble((Double.doubleToLongBits(x) << 1) >>> 1);
    }

    private static double estimate_root(int n, double k) {
        // Extract the exponent from k.
        long exp = (Double.doubleToLongBits(k) & 0x7ff0000000000000L);
        // Format the exponent properly.
        int D = (int) ((exp >> 52) - 1023);
        // Calculate and return 2^(D/n).
        return Double.longBitsToDouble((D / n + 1023L) << 52);
    }
}

Answer 1

Just iterate until the update is less than say, 0.0001, if you want a precision of 4 decimals.

That is, set your epsilon to Math.pow(10, -n) if you want n digits of precision.

Answer 2

Let's recall what the error analysis of Newton's method says. Basically, it gives us an error for the nth iteration as a function of the error of the n-1 th iteration.

So, how can we tell if the error is less than k? We can't, unless we know the error at e(0). And if we knew the error at e(0), we would just use that to find the correct answer.

What you can do is say "e(0) <= m". You can then find n such that e(n) <= k for your desired k. However, this requires knowing the maximal value of f'' in your radius, which is (in general) just as hard a problem as finding the x intercept.

What you're checking is if the error changes by less than k, which is a perfectly acceptable way to do it. But it's not checking if the error is less than k. As Axel and others have noted, there are many other root-approximation algorithms, some of which will yield easier error analysis, and if you really want this, you should use one of those.

Answer 3

You have a bug in your code. Your pow() method's last line should read
return (y & 1) == 1 ? k : k * x;
rather than
return (y & 1) == 0 ? k : k * x;