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