[Bug 277783] libc fma() doesn't not return the correct zero sign

https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=277783

--- Comment #14 from Steve Kargl <kargl@FreeBSD.org> ---
Created attachment 251489
  --> https://bugs.freebsd.org/bugzilla/attachment.cgi?id=251489&action=edit
New patch

The code for fma() and fmal() contain a special case when the addends sum to
zero, but it seems that it mishandles small low order bits especially if a
subnormal condition occurs.  One of the testcases in msun/tests (that I broke
with the original patch) is fma(-1.,1.,1.) with the intermediate results

(x,y,z): -0x1p+0  0x1p+0  0x1p+0
xy = {-0x1p-2, 0x0p+0}
r  = {0x0p+0, 0x0p+0}
spread = 2
zs = 0x1p-2

fma returns (xy.hi + zs +ldexp(xy.lo,spread)), which is fine.

For the case that Victor submitted one has

(x,y,z):  0x1.ffffffffffff8p-501  0x1.0000000000004p-500 -0x1p-1000
xy = {5.000000e-01, -3.944305e-31}
r  = {0.000000e+00 0.000000e+00}
spread = -999 
zs = -5.000000e-01

For the original code, the special leads to
ldexp(xy.lo,spread)=ldexp(-3.944305e-31,-999), is not good.

The new patch skips the special case if xy.lo != 0.  This then allows fma() to
call its add_and_denormalize() function.

-- 
You are receiving this mail because:
You are the assignee for the bug.