Fixing ilogb()

[Bruce Evans]

On Sat, 8 May 2004, David Schultz wrote:

> On Sat, May 08, 2004, Stefan Farfeleder wrote:
> > I found two problems with our current ilogb() implemenation(s).
> >
> > - Both the man page and FP_ILOGB0 claim that ilogb(0) returns INT_MIN
> >   while the implementation in s_ilogb.c returns -INT_MAX.  We really
> >   have to decide on one value (both are conforming to C99, I've attached
> >   the text).  The attached patch assumes that -INT_MAX is kept.
>
> FWIW, SunOS has historically used the following mapping:
>
> 	ilogb(0)        ==> -INT_MAX
> 	ilogb(NAN)      ==> INT_MAX
> 	ilogb(INFINITY) ==> INT_MAX
>
> This matches OS X, our MI ilogb() implementation[1], and your
> patch, and I think those are pretty good reasons to use your
> version.

I agree.

> >   On a related note, is there a reason why <math.h> doesn't use
> >   <machine/_limit.h>'s __INT_M{AX,IN} for the FP_ILOGB* macros?
>
> Yes, machine/_limits.h did not exist when it was written, so there
> was no way to get INT_{MIN,MAX} without namespace pollution.
> It would be a good idea to use these in math.h and s_ilogb[f].c now.

It would also be good to keep using fdlibm's code, but the hard-coded
INT_{MIN,MAX} are too vax-dependent.

> > - On i386 the assembler file s_ilogb.S is used instead.  It uses the
> >   instruction fxtract which returns INT_MIN for 0, infinity and NaN.
> >   Except for infinity this conforms too, but doesn't play well with our
> >   MI definitions for FP_ILOGB*.  In the attached patch I've aligned
> >   s_ilogb.S's behaviour with s_ilogb.c, the alternative would be just
> >   fixing infinity and making FP_ILOGB* MD.
>
> Although it would be legitimate to make FP_ILOGB* MD, we have to
> fix the infinity case anyway, so we might as well make the NAN and
> 0 cases MI.
>
> FWIW, I was working on some faster MD implementations of
> fpclassify() and friends a while ago, and getting these corner
> cases to be consistent is a big PITA.  ;-)

It also has unfortunate runtime overheads.

Returning INT_MIN for the i387 version of ilogb(Inf) was from overflow.
fxtract(Inf) gives Inf, but fistpl(Inf) gives INT_MIN since a preposterous
value is the closest approximation to a NaN that can be represented
in an int.

=== Related bugs ===

gcc still gets most conversions from large values to longs and longs
wrong because it doesn't understand this.  E.g.:

%%%
 * constant 1e25 converted to long long = 0x7fffffffffffffff
 * constant 1e25 converted to unsigned long long = 0x7fffffffffffffff
 * variable 1e25 converted to long long = 0x8000000000000000
 * variable 1e25 converted to unsigned long long = 0x1614014800000000
 * constant Inf converted to long long = 0x7fffffffffffffff
 * constant Inf converted to unsigned long long = 0x7fffffffffffffff
 * variable Inf converted to long long = 0x8000000000000000
 * variable Inf converted to unsigned long long = 0
%%%

Here:
- the first conversion is reasonable (LLONG_MAX is a reasonable approx.
  for Inf as a long long)
- the second not so reasonable (should be ULLONG_MAX to match the first)
- the third is reasonable (LLONG_MIN is a reasonable approx. for NaN as
  a long long) but inconsistent with the first.
- the fourth is garbage
- the fifth, sixth and seventh are consistent with the first, second and
  third, respectively
- the eigth is different garbage than the fourth.

Negative large values and Infs give different behaviour (slightly better,
but not the negatives of the above results).

Conversions to longs and unsigned longs give different behaviour (sightly
worse, starting with conversion of Inf to -1L).

The bugs used to be worse.  Typically, gcc would use fistpq for runtime
conversions to avoid certain problems and then not understand what to
do with the overflow indicator 0x8000000000000000 when converting to
a long.  Compile-time conversions can use infinite precision so they
tend not to be bug for bug compatible, but they now are except for the
fourth and eight cases above.

Of course, all these are not incorrect, since the behaviour is undefined.
I prefer a SIGFPE.

=== End Related bugs ===

% Index: src/lib/msun/i387/s_ilogb.S
% ===================================================================
% RCS file: /usr/home/ncvs/src/lib/msun/i387/s_ilogb.S,v
% retrieving revision 1.8
% diff -u -r1.8 s_ilogb.S
% --- src/lib/msun/i387/s_ilogb.S	6 Jun 2000 12:12:36 -0000	1.8
% +++ src/lib/msun/i387/s_ilogb.S	8 May 2004 18:57:27 -0000
% @@ -43,11 +43,27 @@
%  	subl	$4,%esp
%
%  	fldl	8(%ebp)
% +	fxam
% +	fnstsw	%ax
% +	sahf

This is the main runtime overhead.  I think it can mostly be avoided by
checking the return value.  ilogb() can only be INT_MIN after overflow
or other conversion errors (check this).  There 3 cases:
- logb(0) = -Inf; fistpl(-Inf) = INT_MIN + IE
- logb(Inf) = Inf; fistpl(-Inf) = INT_MIN + IE
- logb(NaN) = same NaN; fistpl(NaN) = INT_MIN + IE
After finding one of these rare cases, the exact case still needs to be
determined by looking at the original value or the result of fxtract.
Then fucom with 0 should be a faster way to do the classification.  A
full classification is not needed sice denormals are not special here
and unsupported formats are unsupported here too.

% +	jc .L2
% +	jnp .L3
% +

Style bugs: space instead of tab.

%  	fxtract
%  	fstp	%st
%
%  	fistpl	-4(%ebp)
%  	movl	-4(%ebp),%eax
% -
% +.L1:

Style bug: the blank line is more needed than before.

%  	leave
%  	ret
% +
% +.L2:
% +	/* Depending on ilogb(NAN) == ilogb(INFINITY) */

Better wordingi and punctuation: "This depends on ... .".  I would
prefer an unambigous value for NaNs.

% +	movl	$0x7fffffff,%eax	/* FP_ILOGBNAN, INT_MAX */

Style bugs: some comments could be improved, but won't be needed when
<machine/_limits.h> is used.

% +	fstp	%st

Style bug: %st(0) is normal style in this directory except in this file
and s_logb,S,

% +	jmp .L1

Style bug: a blank line is needed more here than elsewhere.

% +.L3:
% +	movl	$0x80000001,%eax	/* FP_ILOGB0 */
% +	fstp	%st
% +	jmp .L1

Bruce