Writing Bulletproof Numerical Programs

Can anybody suggest a website or book that thoroughly covers the topic of writing computer programs to properly deal with computer machine error?

All the books I have on the topic of numerical programming talk about the speed of an algorithm’s convergence, the various algorithms available for a particular application, etc. However, what I would like is material that really gets into the bits and bytes of bulletproof numerical programming.

For example, I would like to write a program to compute the roots of a cubic polynomial equation:

ax^3 + bx^2 + cx + d = 0

A technique for finding the exact roots of this equation is available, but I want to ensure that the program is extremely robust. As a first step, I would like to check if a is zero before dividing the entire equation by a:

x^3 + (b/a)x^2 + (c/a)x + (d/a) = 0

Am I correct in simply doing an if statement?



Should I instead be checking





Or should I be checking against yet some other quantity?

I then proceed to compute various quantities along the way to computing the roots of the original cubic equation. Along the way, what kinds of checks should I be making? Against what quantities should I be comparing my computed values? etc.

These are the sorts of questions I have, but cannot find a good resource to answer them.

I will be going through my program line-by-line to make it as bulletproof as possible. Any help in this goal would be much appreciated.

I don't exactly have pearls of wisdom to offer to you, but I'll link to this paper which has achieved status of canonical work. I attended a talk on it given by this guy, a designer of the Java floating-point library.

There was this book from years ago called Numerical Recipes in C that talked in detail about this. See if you can get your hands on a copy.

Neither of your check for zero solutions are good.

A result that is mathematically zero may not be zero in the machine representation because of floating point errors in intermediate calculations.

DBL_MIN is a large magnitude negative number, not zero which is in the centre of the range of double.

DBL_EPSILON is a small positive value. The smallest value that may be added to 1.0 to result in a different stored value. Because it is positive, if anything your conditional should be:Although in practice I would use a value determined by the necessary precision of the algorithim.

You may not need to worry about it. I just tested these: and neither produce a divide by zero exception, but rather 1.#INF00 (which ia a form of not-a-number or NaN).

So a better approach may be to perform the desired operation and test the result using the _isnan( double ) function. This function is not an ISO standard function but an IEEE recommendation. You will probably find it in float.h or math.h depending on your compiler's library implementation.


Clifford

Thanks for the feedback everyone.

Any other resources to suggest?

I have continued to look around and have been browsing the NETLIB website, which is supposed to be a library of excellent-quality numerical programs. Some of the posted programs use quantities similar to DBL_MIN and DBL_MAX but they use them in ways I wouldn’t expect.

For example, say you want to compute the square of a number, x^2. If x is very small, you want to prevent an underflow error. On the other hand, if x is very large, you want to prevent an overflow error. In these cases, I would expect to compare x^2 with DBL_MIN or DBL_MAX directly. Instead, x^2 is compared to (DBL_MIN/DBL_EPSILON), (DBL_MAX/DBL_EPSILON), the square root of DBL_EPSILON, or some other quantity. Any ideas why this seemingly roundabout approach is taken?

Comparing x^2 with DBL_MAX would really do you no good. If x is too big, x^2 will overflow and will still be smaller than DBL_MAX.

In the realm of integers, it's easier to see. Let's assume that we have only 4 bits to work with and let's assume unsigned integers. INT_MIN is 0 (0000) and INT_MAX is 15 (1111). The square of 5 is 25, or 11001. Oops, overflow and we're left with 1001 (9 in decimal). x^2 (9) is smaller than INT_MAX (15), but we know that's not the right answer so that's not a good validity test.

Now, the realm of floating-points is much more complex and I'm not enough of a math-head to explain what tricks are being performed in them there fancy libraries. I would just trust it (blindly ).

Edit: typo.

The problem is actually worse. The C standard (at least the older one) says that an operation resulting in overflow results in undefined behavior. You might know what your particular compiler does, but you can't write portable, robust code that detects an overflow after the event. You need to detect in advance that it's going to happen.

As I thought too. But did you notice that the OP claims this library tests the result (x^2) against certain values? Not only that, but they apparently test against DBL_MAX/DBL_EPSILON, which is an obvious overflow itself. I too certainly would like to see that one explained.

David_B: any chance you can point us directly to the source file. I've browsed the netlib website but it's not straightforward.

Here's a routine from the BLAS library, SNRM2, for computing the Euclidean norm of a real vector.

SNRM2

This routine uses the parameters CUTLO and CUTHI to prevent underflow and overflow, respectively.

CUTLO = maximum of SQRT(U/EPS)
CUTHI = minimum of SQRT(V)

where
EPS = smallest no. such that EPS + 1. .GT. 1.
U = smallest positive no. (underflow limit)
V = largest no. (overflow limit)

I don't know why the square root is taken. In addition, the definitions aren't even the same. Why wouldn't CUTLO simply be SQRT(U) to make its definition analogous to the definition of CUTHI?

In any case, perhaps you could answer one question for me:

When checking a quantity, whether it be user input or computed, against 0, what is the best practice?

Should it use something like CUTLO?

Did you read the linked material regarding floating point?

Time and the fact that it is not my problem we are dealing with means that no I did not. Should I have? Presumably there is something wrong about my statement . You may or may not enlighten me, but right now I am not about to wade through that paper right now. I have bookmarked it though!.


Clifford

See the above comment. Didn't have any points to give .