Legacy Documentclose button

Important: The information in this document is obsolete and should not be used for new development.

Previous Book Contents Book Index Next

Inside Macintosh: PowerPC Numerics / Part 1 - The PowerPC Numerics Environment
Chapter 1 - IEEE Standard Arithmetic


Using IEEE Arithmetic

This section provides some example computations and describes how using IEEE arithmetic in the PowerPC Numerics environment makes programming these computations easier.

Evaluating Continued Fractions

Consider a typical continued fraction cf(x) .

cf(x) = 4+-3 OVER x+-2+-1 OVER x+-7+10 OVER x+-2+-2 OVER x+-3

An algebraically equivalent expression is rf(x) :

rf(x) = 622+-x(751+-x(324+-x(59+-4x))) OVER 112+-x(151+-x(72+-x(14+-x)))

Both expressions represent the same rational function, one whose graph is smooth and unexceptional, as shown in Figure 1-2.

Figure 1-2 Graph of continued fraction functions cf(x) and rf(x)

Although the two functions rf(x) and cf(x) are equal, they are not computationally equivalent. For instance, consider rf(x) at the following values of x:

x = 1
rf(1) = 7
x = 2
rf(2) = 4
x = 3
rf(3) = 8/5
x = 4
rf(4) = 5/2
Whereas rf(x) is perfectly well behaved, those values of x lead to division by zero when computing cf(x) and cause many computers to stop. In IEEE standard arithmetic, division by zero produces an Infinity. Therefore, PowerPC Numerics has no difficulty in computing cf(x) for those values.

On the other hand, simply computing rf(x) instead of cf(x) can also cause problems. If the absolute value of x is so big that x4 overflows the chosen data format, then cf(x) approaches cf( ) = 4 but computing rf(x) encounters (overflow)/(underflow) , which yields something else. PowerPC Numerics returns NaN for such cases; some other machines return (maximumvalue)/(maximumvalue) = 1 . Also, at arguments x between 1.6 and 2.4, the formula rf(x) suffers from roundoff error much more than cf(x) does. For those reasons, computing cf(x) is preferable to computing rf(x) if division by zero works the way it does in PowerPC Numerics, that is, if it produces Infinity instead of stopping computation.

In general, division by zero is an exceptional event not merely because it is rare but because different applications require different consequences. If you are not satisfied with the consequences supplied by the default PowerPC Numerics environment, you can choose other consequences by making the program test for NaNs and Infinities (or for the flags that signal their creation).

Rather than sprinkle tests throughout the program in an attempt to keep exceptions from occurring, you might prefer to put one or two tests near the end of the code to detect the (rare) occurrence of an exception and modify the results appropriately. That is more economical than testing every divisor for zero (since zero divisors are rare).

Computing the Area of a Triangle

Here is a familiar and straightforward task that fails when subtraction is aberrant: Compute the area A(x,y,z) of a triangle given the lengths x,y,z of its sides. The formula given here performs this calculation almost as accurately as its individual floating-point operations are performed by the computer it runs on, provided the computer does not drop digits prematurely during subtraction. The formula works correctly, and provably so, on a wide range of machines, including all implementations of PowerPC Numerics.

The classical formula, attributed to Heron of Alexandria, is

A(x,y,z) = SQRTs(s+-x)(s+-y)(s+-z)

where s = (x+y+z)/2 .

For needle-shaped triangles, that formula gives incorrect results on computers even when every arithmetic operation is correctly rounded. For example, Table 1-2 shows an extreme case with results rounded to five decimal digits. With the values shown, rounded
(x+(y+z))/2 must give either 100.01 or 100.02. Substituting those values for s in Heron's formula yields either 0.0 or 1.5813 instead of the correct value 1.000025.

Evidently, Heron's formula would be a very bad way for computers to calculate ratios of areas of nearly congruent needle-shaped triangles.
Table 1-2 Area using Heron's formula
  CorrectRounding downwardRounding
upward
 
x100.01100.01100.01 
y99.99599.99599.995 
z0.0250.0250.025 
(x+(y+z))/2100.015100.01100.02 
A1.0000250.00001.5813 

A good procedure, numerically stable on machines that do not truncate prematurely during subtraction (such as machines that use IEEE arithmetic), is the following:

  1. Sort x,y,z so that x y .
  2. Test for z x+-y to see whether the triangle exists.
  3. Compute A by the formula

A = SQRT((x+(y+z))(z+-(x+-y))(x+(y+-z)))/4

WARNING
This formula works correctly only if you do not remove any of the parentheses.
The success of the formula depends upon the following easily proved theorem:

THEOREM If p and q are represented exactly in the same conventional floating-point format, and if 1/2 p/q 2 , then p+-q too is representable exactly in the same format (unless p+-q suffers underflow, something that cannot happen in IEEE arithmetic).

The theorem merely confirms that subtraction is exact when massive cancellation occurs. That is why each factor inside the square root expression is computed correctly to within a unit or two in its last digit kept, and A is not much worse, on computers that subtract the way PowerPC Numerics does. On machines that flush tiny results to zero, this formula for A fails because (p+-q) can underflow.


Previous Book Contents Book Index Next

© Apple Computer, Inc.
13 JUL 1996