Exception Handling I: The Critical Example

[David Hough]

In his proposals for floating-point exception handling, W. Kahan uses the
following example of computing a continued fraction and its derivative.
It's a critical case for exception handling proposals because the desired
exception handling for the 0*inf case changes on each iteration.
Here's the underlying code with no exception handling:

void
continued_fraction(N, a, b, x, pf, pf1)
	int             N;
	double         *a, *b, x, *pf, *pf1;
{
	double          f, f1, d, d1, q;
	int             j;

	/*
	 * Evaluate function f at x by continued fraction a[j] and b[j];
	 * function value to *pf, derivative to *pf1
	 */

/*	Assume
		aj finite
		bj finite and nonzero
		x finite

	CRITICAL STEP is evaluating f1= -(d1/d)*q

	If d == 0 and d1 != 0 on step j, then
	on the next step j-1 you want f1 = b(j-1)/b(j) * d1(j),
	a result determined by taking limits as d->0 in step j.
	Instead by IEEE default, 
	on the next step j-1 you will get f1 = -(inf/inf)*0
	or NaN.

	if d == 0 and d1 == 0 on step j, then 
	on the next step j-1 you want f1 = 0 .
	Instead by IEEE default, 
	on this step j you get f1 = -(0/0)*inf
	or NaN which pollutes all subsequent f1's.
*/

	f1 = 0;
	f = a[N];

	for (j = N - 1; j >= 0; j--) {
		d = x + f;
		q = b[j] / d;
		f = a[j] + q;
		d1 = 1.0 + f1;
		f1 = -(d1 / d) * q;		/* CRITICAL STEP */
	}

	*pf = f;
	*pf1 = f1;
	return;
}

Thus IEEE default results for invalid exceptions
are unhelpful for this particular example, producing NaN's for the derivative
f1.  The key point is that the correct result depends on what went on before
in a complicated way.   The straightforward way to code the needed test in
the absence of language support looks like
 
	for (j = N - 1; j >= 0; j--) {
		f1j2 = f1j1;	/* Save this guy for recomputation -
			won't be defined until j=N-2 but that's the
			earliest it will be needed.. */
		f1j1 = f1;	/* Intermediate save. */
		d = x + f;
		q = b[j] / d;
		f = a[j] + q;
		d1 = 1.0 + f1;
		r = d1 / d;	/* CRITICAL STEP */
		if (r != r) { /* r is NaN from 0/0 or inf/inf */
			if (d1 == 0) {
				/* f1 = (0/0)*inf so return infinity */
				f1 = q;
			} else {
				/* f1 = (inf/inf)*0 so return limit value */
				f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];
			}
		}
		else { /* r is a normal number */
			f1 = -r * q;
		}	
	}


The "normal" case is slowed down by the necessity of retaining f1j2,
the value of f1 from two iterations back, which can't be avoided, 
and by having to test
the result of d1/d and branch conditionally.   All the exception-handling
mechanism to be discussed below is intended to avoid that conditional
branch or reduce its cost.  

On non-IEEE systems on which division by zero was fatal or infinity or NaN could
not be distinguished from finite values, producing robust results for this
computation would be much more expensive.  As they also would be on IEEE systems
which stubbornly optimized away if (r!=r), which could also be expressed
as isnan(r) or unordered(r,r) at somewhat higher cost.   
Note that there is nothing
truly exception or abnormal about the case d == 0; all that's required is that
the value of the input parameter x happen to be the negative of one of the
partial function values f, which is mostly a matter of luck.