The Mathematica Journal
Departments
Download This Issue
Home
Feature Articles
Graphics Gallery
Tricks of the Trade
In and Out
Columns
The Mathematica Programmer
New Products
New Publications
Classifieds
Calendar
News Bulletins
Editor's Pick
Mailbox
Letters
Write Us
About the Journal
Staff and Contributors
Submissions
Subscriptions
Advertising
Back Issues

Machine Arithmetic

Q: Why is
[Graphics:../Images/index_gr_131.gif]
[Graphics:../Images/index_gr_132.gif]
different from
[Graphics:../Images/index_gr_133.gif]
[Graphics:../Images/index_gr_134.gif]
A: David Withoff (withoff@wolfram.com) answers: I think that there are two distinct questions here.
1. Why is [Graphics:../Images/index_gr_135.gif] different than [Graphics:../Images/index_gr_136.gif]?
2. Since the exact result is [Graphics:../Images/index_gr_137.gif], how can Floor return 468?
The answer to the first question is that this is a consequence of the limitations of machine arithmetic. The same "bug" is present in nearly all computer software. As distressing as this may seem, the fact is that multiplication of machine numbers is not associative commutative.
[Graphics:../Images/index_gr_138.gif]
[Graphics:../Images/index_gr_139.gif]
[Graphics:../Images/index_gr_140.gif]
[Graphics:../Images/index_gr_141.gif]
The default output form rounds values to six digits, so the numbers appear to be the same. However,
[Graphics:../Images/index_gr_142.gif]
[Graphics:../Images/index_gr_143.gif]
The difference is nonzero and is actually [Graphics:../Images/index_gr_144.gif], a consequence of the way floating point numbers are represented (look at the NumericalMath`Microscope` package, for an example). You can confirm that there is only a one-bit difference between [Graphics:../Images/index_gr_145.gif] and [Graphics:../Images/index_gr_146.gif] using RealDigits.
[Graphics:../Images/index_gr_147.gif]
[Graphics:../Images/index_gr_148.gif]
[Graphics:../Images/index_gr_149.gif]
[Graphics:../Images/index_gr_150.gif]
This difference, although very small, makes a big difference to Floor. Functions like Floor, Ceiling, and Round have discontinuities and so are very sensitive to small perturbations near the discontinuities.
To demonstrate that the same type of behavior exists in other computer software, here is a simple C program


#include <math.h>

main() {
        double x, y;

        x = floor(0.7 * 67 * 10);
        y = floor(67 * 10 * 0.7);

        printf ("floor(0.7 * 67 * 10) = %f\n", x);

        printf ("floor(67 * 10 * 0.7) = %f\n", y);
}

which, on my computer, gives an equivalent result.


floor(0.7 * 67 * 10) = 469.000000


floor(67 * 10 * 0.7) = 468.000000

Lots of people have devoted decades of their lives to stabilizing numerical algorithms against the effects of this and other idiosyncrasies of machine arithmetic.
The answer to the second question, again both in Mathematica and in other systems, is that this is a design decision. An inexact number (a number in which the last digits are uncertain) is usually best thought of as an interval rather than as a point. For example, a number like 2.00000 can be thought of as the best available six-digit representation for any number between 1.9999995 and 2.000005, i.e., Interval[{1.9999995, 2.000005}]. Since this interval includes 2, there is an obvious question about whether the "floor" of that number should be 2 or 1. There are lots of ways to handle this ambiguity and they have different advantages in different situations, but eventually it comes down to a design decision.
In Mathematica, one of the primary reasons for including variable-precision arithmetic is to provide a way to get around the idiosyncrasies of machine arithmetic. Variable-precision is not without idiosyncrasies either, but it is far more reliable and consistent than machine arithmetic. With variable-precision arithmetic rather than machine arithmetic, you should be able to get consistent results from these calculations.
However, the best solution for this and almost any similar problem in numerical analysis, whether in C, Fortran, Mathematica, or any other system, is to arrange your calculation so that it is not sensitive to how the system chooses to handle these ambiguities.


Converted by Mathematica      September 29, 1999

[Prev Page]