Avoiding counting mode with IEEE 754 appendix functions

[Bob Alverson]

> 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.
> 

Sounds good to me.  This operation could be synthesized out of logb
and scalb, if they were fast.  It's not clear that too much effort
should go into handling these cases at the expense of making add and
multiply go fast.  TANSTAFL

> 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.
>

Zero is a little odd for the exceptional cases.  The actual value is
not critical, but they should be outside the range for non-zero numbers,
shouldn't they?  Otherwise one might take logb of a number and conclude
the number is between 1 and 2.

Producing a float result is somewhat artificial, but I can hear the
hardware guys groaning about CPU-FPU transfers from here.  I would
suggest that without that little bit of implementation influencing the
architecture, an integer result from logb would work fine.

> z = scalb(x,y)
> 
> z is x * 2**y computed by exponent manipulation; if y is not an integral
> value, 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.
> 

I would be surprised to see any aggressive FPU implement scalb with a
float power.  The processing is too complex.  Assuming your FPU can
convert float to int and leave the result in a float register, you
could reasonably implement scalb(x, y) as scalb(x, int(y)).  The only
loss is that invalid op exceptions might get turned into
overflow/underflow.  I'm not sure how infinities get into the scaling
factor anyway.

Bob