Avoiding counting mode with IEEE 754 appendix functions

[David Hough]

After looking over several examples (determinant, hypot, Clebsch-Gordan
coefficients) I've come to the conclusion that the most useful primitives
for avoiding gratuitous overflows and underflows are very similar to
those specified in the IEEE 754 appendix:
 
z = significand(x)
 
        If x is normal or subnormal, z is its value scaled to between
1 and B.  If x is zero, infinity, or NaN, z is x unchanged.  significand
raises no exceptions.
 


z = logb(x)

        If x is normal or subnormal, z is the unbiased exponent
(in floating-point format) such that x == scalb(significand(x),logb(x)).
For such x, z is an integral value.
If x is zero, infinity, or NaN, z is zero.  logb raises no exceptions.

z = scalb(x,y)

        z is x * 2**y computed by exponent manipulation; if y is not an integralvalue, then its fraction is discarded (so it is rounded toward zero).
scalb can raise overflow or underflow,
with their associated inexact exceptions, as well as invalid operations such
as 0 * 2**+inf or inf * 2**-inf or use of signaling NaN.


significand and logb raise no exceptions because they are never used in
a vacuum; something else will happen later if it needs to.  scalb is the
thing that happens later, so it must handle exceptions normally.

Using these the inner loop of a determinant calculation is transformed from

        d = d * a(i,i)

to

        ds = ds * significand(a(i,i))
        ld = ld + logb(a(i,i))

and the latter loop is followed by

        d = scalb(ds, ld)


The loop costs four fpops rather than one, but this may be cheaper than
having two conditional branches as in

        d = d * [exception underflow goto label1][exception overflow goto label2] a(i,i)

or an external function call

        d = d * frexp(a(i,i), &la)
        ld = ld + la

and restartable traps, continuable traps, and counting mode are all unneeded.

Another approach is to define some new operators *signif and *logb (and
/signif and /logb) which compute significands and logb's of products and
quotients without actually forming the products or quotients, which might
overflow or underflow:

        d = d *signif a(i,i)
        ld = ld + d *logb a(i,i)

requires three fpops rather than one.  However looking at computing
hypot without gratuitous intermediate overflow/underflow suggests that the
simpler primitives are more versatile.

While logb() "naturally" returns an integer result and scalb "naturally"
accepts an integer second argument, in typical RISC architectures with
separate floating-point and integer data registers it looks like a good
idea to keep everything in floating-point formats... or else add integer
operations on floating-point data registers.

Anyway the hardware question is
how much it would cost to implement significand, logb, and scalb
as hardware (not microcoded) fpops.  The only language question is adding three
more generic intrinsic operators in functional notation.

As for Clebsch-Gordan coefficients (I spelled the name wrong last time),
judging by the formulas on page 1006 of Abramowitz and Stegun, the best
way to compute those would be to have tables of significand(n!)
and logb(n!) and then multiply the significands together, and add the
exponents, using ordinary arithmetic.
Same for coefficients of hypergeometric functions to the extent that
gamma functions can be replaced by factorials.  I'd appreciate comments
from anybody who has done these kinds of calculations.