[PATCH 2/5] math: Optimize frexpf (binary32) with fast path for normal numbers

[Wilco Dijkstra]

Hi Osama,

> +++ b/sysdeps/ieee754/flt-32/s_frexpf.c
> @@ -18,27 +18,36 @@ static char rcsid[] = "$NetBSD: s_frexpf.c,v 1.5 1995/05/10 20:47:26 jtc Exp $";

1995!!! I'd suggest to remove this and the old copyright since your version is
already a significant rewrite (and the final version will have nothing in common).

>+	if (__glibc_likely ((ex - 1) < 0xfe))
>+	  {
>+	    *eptr = ex - EXPONENT_BIAS + 1;
>+	    return asfloat ((hx & (MANTISSA_MASK | SIGN_MASK)) | FREXPF_EXP);
>+	  }

Looks good - but it's worth checking whether return asfloat (hx - *eptr << 23) is
better on most targets (since even 32-bit masks might need multiple instructions).
You can use godbolt.org to compile small samples for many different targets at once.

>      *eptr = 0;
>+	uint32_t ix = hx & 0x7fffffff;
>+	if (ix >= 0x7f800000 || (ix == 0)) return x + x;	/* 0,inf,nan */

This could be written as (ix - 1) >= 0x800000 or even simpler: (int32_t)(hx << 1) <= 0.
The code would be clearer if the *eptr = 0 is done only when we return. And given
zero is more likely than denormals, adding __glibc_likely is a good idea too.

Trailing comments are hard to maintain - I'd remove the one above and update the
earlier comment to only apply to this if statement.

> +	/* Subnormal */

/* Subnormal.  */

> +	x *= two25;
> +	hx = asuint (x);
> +	*eptr = ((hx & 0x7fffffff) >> MANTISSA_WIDTH) - EXPONENT_BIAS - 25 + 1;
> +	return asfloat ((hx & (MANTISSA_MASK | SIGN_MASK)) | FREXPF_EXP);

Try using stdc_leading_zeros (hx << 8) - this gives the correct shift left you need to
normalize hx. It should be faster than the old multiply method for normalization.
To enable benchmarks for this, add separate entries for zero and denormal like:

## name: zero
...
## name: denormal
...

Cheers,
Wilco