Unexpected NANs in complex arithmetic

[Oscar Benjamin]

On 22 June 2016 at 08:14, Antoon Pardon <antoon.pardon at rece.vub.ac.be> wrote:
> Op 22-06-16 om 04:48 schreef Steven D'Aprano:
>> I'm doing some arithmetic on complex numbers involving INFs, and getting
>> unexpected NANs.
>>
>> py> INF = float('inf')
>> py> z = INF + 3j
>> py> z
>> (inf+3j)
>> py> -z
>> (-inf-3j)
>>
>> So far, nothing unexpected has occurred. But:
>>
>> py> -1*z  # should be the same as -z
>> (-inf+nanj)
>
> Also the multiplication of a+bj with c+dj is (ac-bd)+(ad+bc)j
> With your "numbers" this gives.
>
>   (inf*(-1) - 3*0) + (inf*0 + 3*(-1))j
>
> Again the nan part makes perfect sense.

This explains presumably how complex multiplication is implemented
although it's not really a mathematical justification.

The problem is that mathematically complex numbers distribute over
addition. Complex floats (i.e. Python's complex type) for the most
part distribute approximately over addition so that the expansion
Antoon refers to will usually be correct up to floating point error:

   (1)   (x1 + j*y1) * (x2 + j*y2) ~= (x1*x2 - y1*y2) + j*(x1*y2 + x2*y1)

We can think of this as applying the distribution rule twice twice
(once for each pair of brackets). If the first complex number is real
it should have zero real part i.e. y1 = 0. This leads to

   (2)   (x1 + j*y1) * (x2 + j*y2) = x1 * (x2 + j*y2) = x1*x2 + j*x1*y2

which is I think the kind of formula Steve is expecting when being
surprised above. However we also have using the full distributed
formula (1):

    (3)   (x1 + j*y1)*(x2 + j*y2) = x1*x2 - 0*y2 + j*(x1*y2 + x2*0)

If y2 and x2 are real then the RHSs of (3) and (2) are equivalent.
However in your example x2 is infinity.

In essence the problem is that the set of numbers represented by
Python's complex type is the set of complex numbers plus a few
infinities/nans. The special values break the distributivity of
addition so that the algebra implicit in the implementation of the
complex type doesn't hold.

To see that infinities/nans break distribution over addition consider:

    inf = inf*1  = inf*(3-2) = inf*3-inf*2 = inf-inf = nan

If you did want to change these cases to match some particular idea of
the infinities as limits then special casing for pure real or
imaginary complex numbers (e.g. using (2) for complex numbers with
zero imaginary part) would probably lead to the results Steve expects.
I'm not really sure what the gain is though since we can't really do
robust algebra with infinities anyway.

--
Oscar