Differences between complex number implementations in Haskell and Python

Question

I'm trying to map the complex number functionality in Python, to Data.Complex in Haskell, but I've reached a point where they differ, and I am unsure as to why.

In python:

>>> x = 3j
3j
>>> x.real
0.0
>>> x.imag
3.0

In Haskell:

> import Data.Complex
> let j n = 0 :+ n
> let x = j 3.0
> realPart x
0.0
> imagPart x
3.0

So far they look the same. Looks like operating on them doesn't differ much either:

Python:

>>> y = 1 + x
(1+3j)
>>> y.real
1.0
>>> y.imag
3.0

Haskell:

> let y = 1 + x
> realPart y
1.0
> imagPart y
3.0

In isolation + - * / ** all seem to work the same way. However this operation yields two different results:

>>> z = (y - 1) ** 2
(-9+0j)
>>> z.real
-9.0
>>> z.imag
0.0

But in Haskell:

> let z = (y - 1) ** 2
> realPart z
-9.000000000000002
> imagPart z
1.1021821192326181e-15

Why is this?

Answer 1

In Haskell, (**) for Complex is essentially

a ** b = exp (b * log a)

which has many opportunities for bad rounding errors to creep in. (I don't know enough Python to check what it would do with an analogous log-then-exp expression; the thing I tried complained that it wasn't ready to handle log(3j).) It has a bunch of special cases to thwart rounding errors, but none check for a fully-real integer exponent. You might consider this a bug or infelicity and report it to the folks in charge of the Complex type as another special case worth adding to the implementation of (**).

In the meantime, if you know your exponent is integral, you can use (^) (for positive numbers only) or (^^) instead:

Data.Complex> (0 :+ 3) ^ 2
(-9.0) :+ 0.0

Answer 2

Although the results given by the two languages are different, they aren't very different (as others have indicated in the comments). So you might guess that it's just a matter of slightly different implementations -- and you'd be right.

Daniel Wagner indicates that in Haskell, the ** operator is defined as

a ** b = exp (b * log a)

Haskell does some special casing, but most of the time, the operation relies on the general-purpose definitions of exp and log for complex numbers.

In Python, it's a little different: powers are calculated using a polar representation. This approach involves using a different set of general-purpose functions -- most of them basic trigonometric functions over ordinary floating point numbers -- and does almost no special-casing. It's not clear to me that this approach is better overall, but it does happen to give a more correct answer in the particular case you've chosen.

Here's the core of the implementation:

vabs = hypot(a.real,a.imag);
len = pow(vabs,b.real);
at = atan2(a.imag, a.real);
phase = at*b.real;
if (b.imag != 0.0) {
    len /= exp(at*b.imag);
    phase += b.imag*log(vabs);
}
r.real = len*cos(phase);
r.imag = len*sin(phase);

Here, a is the base and b is the exponent. vabs and at give the polar representation of a, such that

a.real = vabs * cos(at)
a.imag = vabs * sin(at)

And as you can see in the last two lines of code, len and phase give the corresponding polar representation of the result, r.

When b is real, the if block isn't executed, and this simplifies to De Moivre's formula. I can't find a canonical formula covering the complex or imaginary cases, but it appears to be pretty simple!