[Jack Jansen]
> It turns out that even simple things like 0j/2 return -0.0.
>
> The culprit appears to be the statement
> r.imag = (a.imag - a.real*ratio) / denom;
> in c_quot(), line 108.
>
> The inner part is translated into a PPC multiply-subtract instruction
> fnmsub fp0, fp1, fp31, fp0
> Or, in other words, this computes "0.0 - (2.0 * 0.0)". The result
> of this is apparently -0.0. This sounds reasonable to me, or is
> this against IEEE754 rules (or C99 rules?).
I've said it twice, but I'll say it once more <wink>: under 754 rules,
(+0) - (+0)
must return +0 in all rounding modes except for (the exceedingly unlikely, as
it's not the default) to-minus-infinity rounding mode. The latter case is
the only case in which it should return -0. Under the default
to-nearest/even rounding mode, and under the to-plus-infinity and to-0
rounding modes, +0 is the required result.
However, we don't know whether a.imag is +0 or -0 on your box; it *should* be
+0. If it were -0, then
(-0) - (+0)
should indeed be -0 under default 754 rules. So this still needs to be
traced back. That is, when you say it computes ""0.0 - (2.0 * 0.0)", there
are four *possible* things that could mean, depending on the signs of the
zeroes. As is, I'm afraid we still don't know enough to say whether the -0
result is due to an unexpected -0 as one the inputs.
> If this is all according to 754 rules the one puzzle remaining is
> why other 754 platforms don't see the same thing.
Because the antecedent is wrong: the behavior you're seeing violates 754
rules (unless you've somehow managed to ask for to-minus-infinity rounding,
or you're getting -0 inputs for bogus reasons).
Try this:
print repr(1.0 - 1e-100)
If that doesn't display "1.0", but something starting "0.9999"..., then
you've somehow managed to get to-minus-infinity rounding.
Another thing to try:
print 2+0j
Does that also come out as "2-0j" for you?
What about:
print repr((0j).real), repr((0j).imag)
? (I'm trying to see whether -0 parts somehow get invented out of thin air.)
> Could it be that the combined multiply-subtract skips a rounding
> step that separate multiply and subtract instructions would take? My
> floating point knowledge is pretty basic, so please enlighten me....
I doubt this has anything to do with the fused mul-sub. That operation isn't
defined as such by 754, but it would be a mondo serious hardware bug if it
didn't operate on endcase values the same way as separate mul-then-sub.
OTOH, the new complex division algorithm may generate a fused mul-sub in
places where the old algorithm did not, so I can't rule that out either.
BTW, most compilers for boxes with fused mul-add have a switch to disable
generating the fused instructions. Might want to give that a try (if you
have such a switch, it may mask the symptom but leave the cause unknown).
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4