Re: [PATCH] math: rewrite fma with mostly int arithmetics

[Rich Felker <dalias@libc.org>]

On Sun, Apr 23, 2017 at 01:00:52PM +0200, Szabolcs Nagy wrote:
> * Rich Felker <dalias@libc.org> [2017-04-22 18:24:25 -0400]:
> > A few thoughts, inline below. I'm not entirely opposed to this, if it
> > turns out to be better than the alternatives, but I would like to
> > understand whether it really is...
> > 
> > On Wed, Apr 19, 2017 at 12:41:40AM +0200, Szabolcs Nagy wrote:
> > > the freebsd fma code failed to raise underflow exception in some
> > > cases in nearest rounding mode (affects fmal too) e.g.
> > > 
> > >   fma(-0x1p-1000, 0x1.000001p-74, 0x1p-1022)
> > > 
> > > and the inexact exception may be raised spuriously since the fenv
> > > is not saved/restored around the exact multiplication algorithm
> > > (affects x86 fma too).
> > 
> > Is it difficult to determine when the multiplication part of an fma is
> > exact? If you can determine this quickly, you can just return x*y+z in
> > this special case and avoid all the costly operations. For normal
> > range, I think it's roughly just using ctz to count mantissa bits of x
> > and y, and checking whether the sum is <= 53. Some additional handling
> > for denormals is needed of course.
> 
> it is a bit more difficult than that:
> 
> bits(a) + bits(b) < 54 || (bits(a) + bits(b) == 54 && a*b < 2)
> 
> this is probably possible to handle when i do the int mul.
> 
> however the rounding mode special cases don't get simpler
> and inexact flag still may be raised incorrectly when tail
> bits of x*y beyond 53 bits are eliminated when z is added
> (the result is exact but the dekker algorithm raises inexact).

One thing to note: even if it's not a replacement for the whole
algorithm, this seems like a very useful optimization for a case
that's easy to test. "return x*y+z;" is going to be a lot faster than
anything else you can do. But maybe it's rare to hit cases where the
optimization works; it certainly "should" be rare if people are using
fma for the semantics rather than as a misguided optimization.

> > > depends on the a_clz_64 patch and previous scalbn fix.
> > > 
> > > fmal should be possible to do in a similar way.
> > > 
> > > i expect it to be faster than the previous code on most targets as
> > > the rounding mode is not changed and has less multiplications
> > > (it is faster on x86_64 and i386), the code size is a bit bigger
> > > though.
> > 
> > Kinda surprising on i386 -- I'd expect the 64x64 multiplications to be
> > costly compared to float ones.
> > 
> 
> i implement 64x64 int mul by four 32x32->64 mul,
> i386 has 32x32->64 mul op so that works out well.

Most archs should have a 32x32->64; if not this is an ISA quality
issue and not one I really want to focus on remedying.

> the float code has to do dekker's exact multiplication
> which uses six(!) fp mul (each of which is probably
> slower than an int mul) and a lot more fp add.

OK, this makes your approach make a lot more sense! Thanks for sharing
that info.

> > > 	/* normalize so top 10bits and last bit are 0 */
> > > 	struct num nx, ny, nz;
> > > 	nx = normalize(ix);
> > > 	ny = normalize(iy);
> > > 	nz = normalize(iz);
> > 
> > If the only constraint here is that top 10 bits and last bit are 0, I
> > don't see why clz is even needed. You can meet this constraint for
> > denormals by always multiplying by 2 and using a fixed exponent value.
> 
> yeah that should work, but i also use clz later

Ah, I missed that. Still it might be a worthwhile optimization here; I
think it shaves off a few ops in normalize().

Rich