In the days when I attended high school, Computer Studies was not even offered as a course. As a result, our universities could not make it a prerequisite for their Computer Science programs. The effect was that students with widely varying levels of programming talent arrived at our hallowed halls of learning, which presented a problem for my old Computer Science Department: If they made CompSci too difficult, then those who had never seen a computer would fail, but if they made it too easy, the hackers would get perfect grades, get bored, not go to class, etc.
The university therefore made the rule that you had to pass Mathematics II before taking Computer Science III. The hardest three years for many a hacker was Mathematics II. We learned a whole lot of things in that course, much of which I never fully understood or appreciated. I found figuring out computer programs much more interesting than reading mathematical proofs.
I wish I had a copy of Dr. Didier Besset's (doctorate in Physics from University of Geneva) book Object-Oriented Implementation of Numerical Methods in those days! It marries numerical methods and programming in a very interesting way. Just look at these algorithms that Dr. Besset covers, implemented in Java and Smalltalk:
- Interpolations: Lagrange, Newton, Neville, Bulirsch-Stoer, and cubic spline
- Zero of function: bisection algorithm, Newton's method, roots of polynomials
- Integration of functions: trapeze, Simpson, Romberg
- Series: infinite series, continued fractions, incomplete gamma function, incomplete beta function
- Linear algebra: vectors, matrices, and all that you might dream of for linear algebra
- Elements of statistics: moments, histograms, random number generators, probability distributions
- Statistical analysis: Fisher-Snedecor, many others
- Optimization: hill-climbing, Powell's algorithm, genetic algorithm
- Data mining: data server, covariance, Mahalanobis distance, cluster analysis
It cuts, it slices, it dices! If seeing all those algorithms doesn't excite you, then you're in the wrong profession. ;-)
Dr. Besset sets the scene in the book by pulling out some performance statistics. We all know Java is slow, right? Have a look at these stats from the book:
| Operation | Units | C | Smalltalk | Java |
| Polynomial 10th degree | msec. | 1.1 | 27.7 | 9.0 |
| Neville Interpolation (20 points) | msec. | 0.9 | 11.0 | 0.8 |
| LUP matrix inversion (100 x 100) | sec. | 3.9 | 22.9 | 1.0 |
The C measurements were done using published algorithms—Dr. Besset didn't just add a whole lot of wait statements into the C code. He writes, "I want to emphasize here that all the code in this book is real code that I have used personally in real applications." Wow! That's certainly better than most books nowadays. :-)
Comparing Doubles Explained
Besides all the interesting algorithms, which are shown with mathematical explanations and well-written Java code, Dr. Besset also tackles issues such as the problems associated with comparing floating-point numbers. Here's a sneak peek at part of an article he recently wrote for the Smalltalk Chronicles (edited by me):
|
One classical caveat with floating-point numbers is checking the equality between two floating-point numbers. Now and then one bloke complains on some news groups that Smalltalk does not compute right with floating-point numbers. In the end it turns out that he was computing a result with method A, the same result with method B, and, to check the results, was evaluating the expression ' resultA == resultB. The fact that this expression evaluates to false has nothing to do with Smalltalk. It is a fundamental problem with floating-point numbers [HK: also in Java].
A floating-point number is only an approximation of a mathematical real number. A small introductory article like this one is too short to explain things in depth, but I would like to quickly recall a few principles.
Floating-point numbers are used to keep the relative error constant. This is valid of course for a given number. As soon as numbers are combined together one must follow the propagation of rounding errors. Because the relative error is kept constant, nothing serious happens with multiplications and divisions. The error on additions and subtractions, however, can become prohibitively large, up to the point of generating something utterly wrong. To illustrate this point, try running the following Java program:
public class DoubleTest {
public static void main(String[] args) {
System.out.println(2.71828182845905 - 2.71828182845904);
System.out.println(2.71828182845905 - 2.71828182845904 +
0.00000000000001 );
System.out.println(2.71828182845905 == (2.71828182845904 +
0.00000000000001));
}
}
[HK: The following is the answer]
9.769962616701378E-15
1.976996261670138E-14
false
Notice the difference in the last digits. The result you will get is 100% wrong!
Unless you are a very good and courageous mathematician, I would not recommend you to attempt to predict error propagation. The easiest and surest thing to do is to measure error propagation experimentally.
After coding an algorithm, you can predict roughly where the infinities or the nearly zero cases are located. I am not speaking only about the result. All steps of the algorithm must be checked against the occurrences of infinities. In these areas, try to evaluate a few results by changing the values by a very small amount (10-12 or so for standard IEEE double format). In general, the difference between the results will be one or two order of magnitudes larger than the original variation. If you observe something much larger, the algorithm used is not made for computers and must be adapted. I give several examples of such modifications in my book. |
How Should You Compare Doubles?
Over a year ago, a friend mentioned that he was surprised that Java programmers did not know how to compare doubles. If you just use "==" as in our example above, you will get incorrect results. Dr. Besset also has a section about that in his book. It is my understanding that in Java the precision of doubles and floats is defined by the IEEE 754 floating-point format, so physical architectures should not have any differences between them, since in Java we are running on a virtual machine. The results were identical on Wintel, an AIX box running Java 1.3 on AIX version 4.3.3.0 (IBM 2 processor 4 cpu model 7044-270), and on a Solaris Box (SunOS Mars2 5.6 Generic_105181-29 sun4u sparc SUNW,Ultra-5_10). The fact that the results are the same can lull you into a false sense of security—they are not guaranteed to be the same, unless you use the strictfp keyword, but then you could lose on performance.
Dr. Besset presents the following algorithm for comparing double-precision numbers (reproduced with permission):
/**
* This class determines the parameters of the floating point
* representation
* HK: I have refactored it somewhat to make it thread-safe and
* to make it easier to understand and to fit into my newsletter.
* The algorithms are the same as in the book.
*
* @author Didier H. Besset
*/
public final class DhbMath {
/** Radix used by floating-point numbers. */
private final static int radix = computeRadix();
/** Largest positive value which, when added to 1.0, yields 0 */
private final static double machinePrecision =
computeMachinePrecision();
/** Typical meaningful precision for numerical calculations. */
private final static double defaultNumericalPrecision =
Math.sqrt(machinePrecision);
private static int computeRadix() {
int radix = 0;
double a = 1.0d;
double tmp1, tmp2;
do {
a += a;
tmp1 = a + 1.0d;
tmp2 = tmp1 - a;
} while (tmp2 - 1.0d != 0.0d);
double b = 1.0d;
while (radix == 0) {
b += b;
tmp1 = a + b;
radix = (int)(tmp1 - a);
}
return radix;
}
private static double computeMachinePrecision() {
double floatingRadix = getRadix();
double inverseRadix = 1.0d / floatingRadix;
double machinePrecision = 1.0d;
double tmp = 1.0d + machinePrecision;
while (tmp - 1.0d != 0.0d) {
machinePrecision *= inverseRadix;
tmp = 1.0d + machinePrecision;
}
return machinePrecision;
}
public static int getRadix() {
return radix;
}
public static double getMachinePrecision() {
return machinePrecision;
}
public static double defaultNumericalPrecision() {
return defaultNumericalPrecision;
}
/**
* @return true if the difference between a and b is less than
* the default numerical precision
*/
public static boolean equals(double a, double b) {
return equals(a, b, defaultNumericalPrecision());
}
/**
* @return true if the relative difference between a and b is
* less than precision
*/
public static boolean equals(double a, double b, double precision) {
double norm = Math.max(Math.abs(a), Math.abs(b));
return norm < precision || Math.abs(a - b) < precision * norm;
}
}
Object-Oriented Implementation of Numerical Methods has details as to why the algorithms work the way they do. Here is how you would use the equals() method:
public class BetterDoubleTest {
public static void main(String[] args) {
System.out.println("Floating-point machine parameters");
System.out.println("---------------------------------");
System.out.println("radix = " + DhbMath.getRadix());
System.out.println("Machine precision = " +
DhbMath.getMachinePrecision());
System.out.println("Default numerical precision = " +
DhbMath.defaultNumericalPrecision());
System.out.println(DhbMath.equals(
2.71828182845905,
(2.71828182845904 + 0.00000000000001)));
System.out.println(DhbMath.equals(
2.71828182845905, 2.71828182845904));
System.out.println(DhbMath.equals(
2.718281828454, 2.718281828455));
System.out.println(DhbMath.equals(
2.7182814, 2.7182815));
}
}
On the machines that I ran this test on, the output was:
Floating-point machine parameters
---------------------------------
radix = 2
Machine precision = 1.1102230246251565E-16
Default numerical precision = 1.0536712127723509E-8
true
true
true
false
I hope you will consider these issues when next you want to compare doubles. For more information on these and other interesting topics, I would recommend that you get a hold of Didier Besset's book.
 | Dr. Heinz Kabutz owns Maximum Solutions, a Cape Town-based consulting firm that specializes in Java technology. He spends the majority of his time programming Java OO applications and also advises companies who wish to embrace Java as a technology. He has been a lead programmer in one of the first big Java developments in South Africa, now consisting of almost 600,000 lines of Java code. He has been programming in Java since 1997. He can be reached at h.kabutz@computer.org. |
