Condition of square root algorithm (Newton's method)

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}\$$

• The first line of your code looks very, very wrong. – gnasher729 Oct 5 '18 at 12:08
• Well, it teaches students not to trust a textbook without question, so that is a good thing :) I hope they also show a good implementation in the book, because this one seems awful. – gnasher729 Oct 5 '18 at 12:22
• I just notice that for denormalized c, this is likely an infinite loop. Shame on Sedgewick and Wayne. – gnasher729 Oct 5 '18 at 12:29
• @Apass.Jack, thanks for clarifying which textbook it is. Curiously this code doesn't seem to be included in the official GitHub repo. – Peter Taylor Oct 5 '18 at 12:51
• @gnasher729 Oh sorry, that's a mistake. I have a eBook and a paper book. The code from the paper book is right. – Rick Oct 5 '18 at 13:02

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.)

• Yes, I should use relative value as you describes. I tried 1e-30 and 100000 with 1e-16, it works fine. So I don't understand why < 1e-15 would cause error. – Rick Oct 6 '18 at 5:40
• Btw, I tried with sqrt(100) and I get 10 exactly, I mean, no decimal. That kind of surprised me. See screenshot here imgur.com/2vAUb8Q. Wow, t + c/t equals 20 exactly. – Rick Oct 6 '18 at 5:44

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.

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).

• I'm not sure how this answers the question. What would go wrong if we used this stopping condition? Also, what is ULP? – Yuval Filmus Oct 5 '18 at 7:28
• @YuvalFilmus Ya.. I also go check what ULP is. I thought it might be a math abbr and then found it seems to be a Java function from wikipedia. But I don't find anything useful there. wikiwand.com/en/Unit_in_the_last_place – Rick Oct 5 '18 at 7:33
• @Rick, I'm not familiar with the textbook you're using, but if it teaches numerical algorithms without teaching such elementary concepts of numerical algorithms as ULPs you're better off finding a different textbook. – Peter Taylor Oct 5 '18 at 7:43
• @PeterTaylor I am reading Algorithm 4th and the Newton formula comes from my mathematical analysis book. Wow, that would be a burden if I have to review all my math books before understanding this.. I mean, I just want to get some intuitive answer or opinion here.. Thank you anyway. – Rick Oct 5 '18 at 7:53
• Ulp = value if lowest bit in the mantissa of a floating point number. – gnasher729 Oct 5 '18 at 12:10