Exception Handling III: Operator Modifiers

[David Hough]

In my first posting, I described a loop to evaluate a continued fraction
and its derivative and how it might be made robust in the face of 0/0 or
inf/inf by conventional means costing a conditional branch on every iteration.

In my second posting, I described how presubstitution might work.  In the best
possible implementation, an extra fp multiply, divide, and "presubstitute"
are required on every iteration.   The "presubstitute" might be a load to
a hardware presubstitution register or a store to a user-space memory location.

I now turn to another approach which is also based on conditional branches -
and thus is probably no cheaper than the conventional approach in this
example.   The payoff is in other situations where the exceptionalness
can't be so readily determined by doing just one floating-point comparison,
such as when an inexact must be detected.

The idea is that instead of testing the numerical operands or result of an
operation, or depending on hardware presubstitution registers or recovering
from an asynchronous trap, the code tests the accrued IEEE exception bits.
It's intended to permit an important optimization on IEEE systems that
maintain separate "current exception bits" for the last fpop: if there is
only one fpop that can generate the exception of interest, then only the
current exception bits for that fpop need be tested.  Abstractly this is not
inherently more or less synchronizing than a conditional branch on 
just-computed floating-point data, although current implementations might
do a better job on the latter.

Conventionally, if you want to test an fp expression for a particular IEEE
exception such as FE_INVALID, you would do something like

	tinv = get_accrued(FE_INVALID);
	clear_accrued(FE_INVALID);
	dest = expression;
	tinv2 = get_accrued(FE_INVALID);
	set_accrued(FE_INVALID, tinv | tinv2);
	if (tinv2 == 0) {
		/* code for normal case that no exception arises */
	} else {
		/* code for unusual case that exception arises */
	}

So that's at least four function calls, worst case, resulting best case in
two stores and two loads of the accrued exception register; plus a conditional
branch.   On most current
high-performance IEEE systems those loads, stores, and branch will be much more
expensive than the "dest=expression" that one would most typically want to
test.

The proposal below is for language support for the foregoing construct, with
the intent that it would often be optimized to 

	dest = expression;
	if (get_current(FE_INVALID) == 0) ... else ...

and that consequently there would in time be hardware support for either the
five instructions

	"branch conditionally if the previous fpop generated FE_X"

for each IEEE exception INVALID, OVERFLOW, etc., or else one instruction

	"branch conditionally if the previous fpop generated FE_X|FE_Y|..."

where any subset of the five IEEE exceptions can be OR'd together
to specify the desired condition; and thus after further time fpop
implementations might send early information to branch units that certain
IEEE exceptions are guaranteed not to arise on this fpop code - a MIPSy idea
much appreciated by Kahan - releasing control flow past the conditional 
branch.

In my experience the most common case, by far, is that only one operation
matters, so my original idea was to attach the branch condition as a modifier
directly to the operation:

	z = x / y ;

becoming

	z = x / [on FE_INVALID goto label] y ;

indicating that the fpop(s) implementing / be tested for invalid and when
encountered, control would flow to a local labeled statement.   
A goto was chosen with malice aforethought to insure a lightweight
implementation burden: if the exception arises the compiled code is neither 
required to have assigned to z nor to have left z untouched.
There is no possible way to return back into the middle of the expression.

A semantically similar but syntactically quite different approach, inspired by
Stallman, is to attach a cast-like modifier to an expression:

	(__oninvalid__ label) (x/y)

or maybe even

	(__oninvalid__ label) (z=x/y) 

The distinction between these two
is important on extended-precision-based systems that
typically won't detect overflow or underflow on the divide but only on the
assignment forcing conversion to storage precision.

Anyway if the tested code is an expression rather than a single operation,
compilers should eventually recognize possibilities for optimization.

An even more general approach is inspired by C++:

	try { z=x/y ; }
	catch (FE_INVALID) { /* code for exceptional case */ } ;
	/* followed by code for non-exceptional case 
	   that will also be executed if exceptional case falls through */

There can be multiple catch clauses for different exceptions.

I don't think the C++ try/catch mechanism is appropriate for C - too
heavyweight: it allows for user-defined exceptions and requires dynamic
scoping of handlers: if FE_OVERFLOW is not caught here, is should be caught
by any dynamically enclosing catch even in some other procedure.   This is
too much like a long jump - costly and not helpful for what I want to do,
which is very local.

Worse, the main purported advantage - 
a cleaner syntax devoid of gratuitous gotos and labels -
will frequently prove illusory, since it will frequently be the case
that exceptional-case code should not fall through to normal-case code as
it would appear to do in the example above.

Thus a truly satisfactory syntax awaits discovery and promulgation.
I'm eager for good suggestions.

I'll close by illustrating the foregoing with the continued fraction example
from the previous postings.   First with the operator or expression
modifier syntax:

void
continued_fraction(N, a, b, x, pf, pf1)
	int             N;
	double         *a, *b, x, *pf, *pf1;
{
	double          f, f1, d, d1, q, r, f1j1, f1j2;
	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 = p = b(j-1)/b(j) * d1(j), a result determined by taking limits
	 * in step j. Since f1(j-1) is nominally computed as (inf/inf)*0, the
	 * desired result can be obtained by substituting inf for inf/inf and
	 * p for inf * 0.
	 * 
	 * if d == 0 and d1 == 0 on step j, then on the next step j-1 you want
	 * f1 = p = 0 . Since f1(j) is nominally computed as (0/0)*inf, and
	 * since f1(j-1) is nominally computed as (f1(j)/inf)*0, the desired
	 * result can be obtained by substituting inf for 0/0 and p (= 0) for
	 * inf * 0.
	 */

	f1 = 0;
	f = a[N];

	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;
#ifdef operatorsyntax
		r = d1 / [on FE_INVALID goto invdiv] d;
#else				/* cast expression syntax */
		r = (__oninvalid__ invdiv) (d1 / d);
#endif
		f1 = -r * q;	/* non-exceptional case */

endofloop:	{		/* some kind of null statement required 
				   syntactically here */
		}
	}

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

invdiv: /* invalid division 0/0 or inf/inf detected */

	if (d1 == 0) {
		/* f1 = (0/0)*inf so return infinity */
		f1 = q;
	} else {
		/* f1 = (inf/inf)*0 so return "p" */
		f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];
	}
	goto endofloop;
}


The corresponding try/catch looks a little nicer:

void
continued_fraction(N, a, b, x, pf, pf1)
	int             N;
	double         *a, *b, x, *pf, *pf1;
{
	double          f, f1, d, d1, q, r, f1j1, f1j2;
	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 = p = b(j-1)/b(j) * d1(j), a result determined by taking limits
	 * in step j. Since f1(j-1) is nominally computed as (inf/inf)*0, the
	 * desired result can be obtained by substituting inf for inf/inf and
	 * p for inf * 0.
	 * 
	 * if d == 0 and d1 == 0 on step j, then on the next step j-1 you want
	 * f1 = p = 0 . Since f1(j) is nominally computed as (0/0)*inf, and
	 * since f1(j-1) is nominally computed as (f1(j)/inf)*0, the desired
	 * result can be obtained by substituting inf for 0/0 and p (= 0) for
	 * inf * 0.
	 */

	f1 = 0;
	f = a[N];

	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;
		try { r = d1/d ; }
		catch (FE_INVALID) { 	/* invalid division 0/0 or inf/inf 
					   detected */
			if (d1 == 0) {	/* f1 = (0/0)*inf so return infinity */
				f1 = q;
			} else {	/* f1 = (inf/inf)*0 so return "p" */
				f1 = (1.0 + f1j2) * b[j + 2] / b[j + 1];
			}
			goto endofloop;
		};

		f1 = -r * q;

endofloop:	{		/* some kind of statement required here */
		}
	}

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

By the way, figuring out that the important exceptional cases are
(0/0)*inf and (inf/inf)*0,  given the input constraints on x, a, and b,
is an interesting exercise.