Condition of square root algorithm (Newton's method)

Question

I am reading the book Algorithm 4th and see this square root algorithm using Newton's method (java sample code):

public static double sqrt(double c)
{
    if (c < 0) return Double.NaN;
    double err = 1e-15;
    double t = c;
    while (Math.abs(t - c/t) > err * t)
        t = (c/t + t) / 2.0;
    return t;
}

I wanna know what's the rule behind the while loop condition,

Math.abs(t - c/t) > err * t

Why is it (|t-c/t|)/t > err? For exmaple, can't the condition simply be |t-c/t| > err?

I know that, t-c/t -> 0 when loop count -> ∞.

I remember that I encouter many "manually choose approximation condition" when learn Probability Theory. But now I totally forget it. So any related,basic background knowledge explanation or introduction would be much appreciated.

---

Newton's method from my original math textbook:

$$ a>0, x_0>0, \\ \{X_n\}: X_n = \frac{1}{2}(x_{n-1} + \frac{a}{x_{n-1}}) \\ \lim{x_n} = \sqrt{a}\ $$

Answer 1

Firstly, there is a typo in the code. if (c > 0) return Double.NaN; should be if (c < 0) return Double.NaN;. To be fair to the questioner, this typo exists in the original textbook, Algorithms, fourth edition by Robert Sedgewick and Kevin Wayne (the e-book version only? I do not have the paper version, yet). To be fair to the book, just one page before where the typo occurs, the same code appears correctly with more explanation on Java syntax.

Secondly, there is a minor magic about the code. Note If you let c=0 and hence t=0 initially, the code c/t will not throw "java.lang.ArithmeticException: / by zero", which you would encounter if running java code int oneOverZero = 1/0;. The call to sqrt(0) will return 0 correctly. This is the magic by IEEE-754 (1985) floating-point standard.

---

Now let us get to the real business.

Math.abs(t - c/t) > err * t

This code is written according to general best practice for checking if two floating-point numbers is equal across most if not all programming languages that implement IEEE Floating-Point Standard, including c/c++ and Java. For some basic background knowledge explanation or introduction, you can checked Comparing Floating Point Numbers, 2012 Edition, especially the section "Relative epsilon comparisons". You can also check the Floating-point Comparison by C++ boost libray.

In a nutshell, in order to find if two floating-point numbers are equal or near enough, instead of comparing them using the "==" operator, or checking whether their difference is 0 or near 0, we should check whether their absolute difference is small enough relative to the minimum (or the maximum) of the two especially when those two numbers are near 0 or very large.

Let us see what might happen otherwise. Suppose we have the following code instead, where we check the absolute difference of t and c/t is smaller than err or not.

public static double sqrtDubious(double c)
{
    if (c > 0) return Double.NaN;
    double err = 1e-15;
    double t = c;
    while (Math.abs(t - c/t) > err)  // is the absolute diff almost 0?
        t = (c/t + t) / 2.0;
    return t;
}

Then sqrtDubious(100000) will run into infinite loop while sqrtDubious(1e-30) returns an apparently wrong number 1.235E-15. However, sqrt(100000) and sqrt(1e-30) works. In fact, as far as I have tried, sqrt(double c) works pretty well for all numbers within the range of Java type double.

Also note the value of err is critical, too. If you change it to be smaller than 1e-15 such as 1e-16, then sqrt(double c) may behave erratically. Sometimes it goes into infinite loop, sometimes not. An accurate description of these wrong cases is rather too hard.

By the way, my testing code is written in Java, the language in which the code in the question is written. So, if you see a different result with your C/C++ code on your machine in the cases where I claim the result is wrong, I will not be too surprised since all these errors cases are very context-sensitive and seemingly erratic, although there are certain patterns if we digger deeper. (Various C/C++ doubles may not be the same as the Java double, anyway.)

Answer 2

Try calculating the square root if 1e-100 with the changed loop condition and see what result you get. Then try calculating the square root of 1e100. Do you run into any problems?

This is numerical mathematics, limits won’t work well. Print out the value of t -c/t and check whether it actually converges, and check if it converges to 0. For some different values if c.

From a wiser man than me: The problem are not the things we don’t know. The problem are the things we know that ain’t so.

Answer 3

can't the condition simply be |t-c/t| > err?

Only if you restrict c such that err will always be more than 1 ULP. Otherwise the condition is equivalent to |t-c/t| > 0, and rounding errors mean that oscillation around the root is possible.

(In actual fact, the chosen value of err * t is about 4.5 ULPs).