more on presubstitution

[David G. Hough at validgh]

I received a message from Dennis Brzezinski that got me thinking about
presubstitution again.    The simplest idea I came up with was:

Have four hardware modes per exception:

	IEEE nonstop default
	nonstandard nonstop default (uflo->0, oflo, div0 -> huge)
	quash store to destination register
	asynchronous non-continuable trap

The first and last you are familiar with.    The second is to provide for
those who regret the passing of the VAX and 370: operations on normal
and zero operands produce only normal and zero results, but it's otherwise
like 754 default.   This mode provides global changed exception handling.

The third is for implementing presubstitution.   The idea is that
you load the destination register with the result to be substituted in
case of exception X, then execute the desired instruction with exception X
in quash store mode.   If exception X arises, the store to the destination
register is quashed.    
This requires a compiler discipline that avoids
using the same register as source and destination.
This is specifically intended to partly 
solve a problem posed by Kahan of computing a continued fraction and its
derivative,
in which the presubstituted result changes at each step of
an iteration:

/*      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.
*/
 
        presubstitute("0/0", infinity());
        presubstitute("inf/inf", infinity());
 
        f1 = 0;
        f = a[N];
        bj1 = b[N];
 
        for (j = N - 1; j >= 0; j--) {
                d = x + f;
                d1 = 1.0 + f1;
                bj=b[j];
                q = bj / d;
                f1 = -(d1 / d) * q;             /* CRITICAL STEP */
                f = a[j] + q;
                p = bj1 * d1 / bj;
                presubstitute("0*inf", p);
                bj1 = bj;
        }

But this example illustrates the full nastiness of the problem.
The invalid exceptions 0/0 and inf/inf need to be handled in a way that
is global, is different from IEEE defaults, and is different from the
invalid exception 0*inf, which must be handled differently on each iteration,
should it arise.    

Anyway, after you've thought of fast economical ways to handle this example,
you can try to think of a way of implementing counting mode to handle
long series of multiplications, e.g. Clebsch-Gordan coefficients,
without getting fouled up by intermediate overflow and underflow.

==========================================================================

Now I remember why my previous musings on this topic
were inconclusive.    For those who have joined this list in the last couple
of years, here are some of them again: