Title: More Notes on Code Critique Competition 12

Francis Glassborow's Student Code Critique Competition is a section of C Vu that I ensure I read in every issue. I can usually make a start at them, to the point of knowing what I would do next to isolate the problems.

But sticking my neck out and attempting to complete them is something I prefer to leave to those with more time and understanding of 'C' so you will generally not find any response from me to this particular column. And anyway, those who wish to take a swing at me will find quite enough material of mine with my neck over extended elsewhere in the pages of this Journal.

So Competition Number 12, on the C++ example finding square roots by Newton's method has been set, the responses printed in the December 2001 issue of C Vu and the prize awarded. Yet I think the evidence is incomplete on one point. So perhaps I will write a few words on this particular competition.

Let me first say that I know a little of 'C', but nothing of 'C++'. The example is in the latter, but I wish to make my comments with respect to the former. So I hope that someone will say "all you have to do to correct this in 'C++' is...". If that someone responds that there is no fix in 'C++', then I think the Standards Committee needs another session.

The problem is the line:

const double TOLERANCE = 1.0e-7;

This variable is used to detect when the difference between the square of the square root estimate and the number being rooted is so small that further calculations will not further refine the estimate. This happens when that difference approaches the accuracy of the floating-point numbers on the target machine. To further refine the estimate would require numbers more accurate than a floating-point number can hold. And a pretty good estimate of that point on the usual 32-bit floating point numbers is, of course, 1.0e-7 which is the value used.

And as The Harpist observes, this value seems pretty generous given that the example uses double precision floating point numbers throughout. A value of 1.0e-15, or perhaps slightly smaller, would be more appropriate in this case.

But who said anything about the example using 32-bit single precision, or 64-bit double precision, numbers?

ICL could process 24 bit floating-point numbers when asked. All right, those were quite limited and would now be prohibited by the 'C' standard. At the other end of the scale, CDC (a long since deceased supplier of behemoths) used 60 bits for their floating point numbers. With Intel working on their IA-64 architecture, we cannot be far off the day when every desktop machine uses 64 bits for single precision and 128 bits for double precision floating point numbers. Some workstations already do. And as for long double...

So how do you find out what a sensible value for the TOLERANCE variable is on your machine?

The 'C' standard (and this is where I switch into 'C', and invite someone to extend this piece to 'C++') has thought of this problem and specified two header files: <limits.h> and <float.h>. They tell you useful things about the numbers on the computer you are using. <limits.h> deals with the sizes of, and largest and smallest numbers that can be held in, character and integer variables. For instance, this is the place to look if you want to know how many bits there are in a byte: it may be more than 8. Look at CHAR_BIT to see what your compiler does on your machine.

But I digress: the file of interest in this problem is <float.h>, which tells us all about how the compiler implements floating point numbers, of whatever length, on the target hardware. And FLT_EPSILON gives the relative accuracy to which numbers can be stored in a single precision floating point number, which is one of type float. DBL_EPSILON and LDBL_EPSILON do the same for double and long double floating point numbers respectively.

The temptation is to initialise TOLERANCE with the value of DBL_EPSILON and then, as The Harpist says, scale this by x which is the number whose square root is to be found.

This is necessary because all the *_EPSILON symbols are actually "the difference between 1.0 and the least number greater then 1.0 that is representable in the given floating point type." Or, if you prefer, the difference between 1.0 and the next largest number. So DBL_EPSILON says the sort of accuracy you will get near 1.0 for arithmetic with double, you have to multiply it by x to find the accuracy you can expect near x.

But to do so may, and probably will, result in the loop in sqroot failing to terminate for certain arguments. Why? Because the computer hardware can only approximate the result of each operation required to form the next estimate of the square root. That is, each operation introduces a small error into the new estimate. So the square of the estimated square root may never get as close to x as x * DBL_EPSILON.

Newton-Raphson iteration, which is the generalisation of Newton's method, is usually well behaved when finding a square root, so this extra error introduced by the computer is normally unimportant. But when the estimate is close to the value of the square root, this can result in successive estimates which oscillate either side of the true value, but neither of which is within x * DBL_EPSILON of the true value.

I spent Sunday morning trying to find an example of this happening and failed. However, if you try the approach of sqroot with integer arithmetic (and remember that on divide, integer arithmetic always discards the remainder so that the result is truncated towards zero) you will get results with the numbers 1022 and 1023 which demonstrate the sort of problems that can arise. Both of these numbers have square roots just less than 32.

Trying to find the square root of 1023 will give estimates which alternate between 31 and 32. So a TOLERANCE of 1, which is the 'epsilon' for integer arithmetic, will stop the iteration with a value of 32. But not 31.

With an argument of 1022, the approach yields a square root of 31, which does not change with further iteration, and whose square is 961. The error is 61 and no sensible value of TOLERANCE will detect that so the loop will run forever.

The solution is to make TOLERANCE a little bigger. If n arithmetic operations are required to form each estimate (and n is 5 here: an add, two divides, a multiply and a subtract), then somewhere between sqrt(n) and n times DBL_EPSILON should be sufficient.

So the original line must be replaced with the two lines of:

#include <float.h>
const double TOLERANCE = 5 * DBL_EPSILON;

I have chosen to multiply by 'n' here.' If you want a more accurate approximation, you could try a smaller multiple at the risk of the iteration never terminating. You can stop that by counting the number of iterations and stopping when this becomes too large. You can try predicting what 'too large' means: I would guess more than 10 iterations for 'double' arithmetic, but you really have to try it and see.

As a finale, if you really want to evaluate a square root, or any other mathematical function on a computer, about the only book which will tell you how to is the Software Manual for the Elementary Functions by William J Cody, Jr, and William Waite, published by Prentice-Hall, 1980, ISBN: 0-13-822064-6.

It reads more like a cookbook than a textbook, and has no mathematics at all in it, just flow charts with careful explanations for you to code. This is a 'how to' book, not a 'why' book. The book tells you how to check that you coded the flow chart correctly for your chosen function by checking both some random general, and specified critical, cases. Interestingly, it does this by presenting a FORTRAN listing, not another flow chart.

My only quibble is that this book does now need another edition to cater for the larger precision computers that are starting to appear which provide floating point numbers of 128 bits or more. Such an edition could also address the "photo-reduced typescript" look of the current edition which was caused by the production techniques in vogue when it was published.

As I now seem to have slipped a book review into the 'Features' section, I am going to stop before the Editor notices. But please keep those Code Critiques coming: I have never found as many points in the problem set as any of your published Code Critiques.