[Python-Dev] test_coercion failing

Tim Peters tim.one@home.com
Wed, 21 Mar 2001 15:18:54 -0500


[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).