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

	Correct	Rounding downward	Rounding upward
x	100.01	100.01	100.01
y	99.995	99.995	99.995
z	0.025	0.025	0.025
(x+(y+z))/2	100.015	100.01	100.02
A	1.000025	0.0000	1.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.