Data.Complex package gives strange result for cube root of -8 + 0i ((-8):+0)

Question

I have to say, math is not my strong suit. I was hoping to get a decent result for How to find the cube root of a negative integer such that it does not return NaN? by using the Data.Complex package but when i do like

*Main> ((-8):+0) ** (1/3)
1.0 :+ 1.732050807568877

I expected to get something like (-2.0):+0 where -2 is the real part and 0 is the imaginary part. However the result turns out with a significant imaginary part. I have checked the (**) RealFloat instance of the Complex type where it states;

x ** y = case (x,y) of
  (_ , (0:+0))  -> 1 :+ 0
  ((0:+0), (exp_re:+_)) -> case compare exp_re 0 of
             GT -> 0 :+ 0
             LT -> inf :+ 0
             EQ -> nan :+ nan
  ((re:+im), (exp_re:+_))
    | (isInfinite re || isInfinite im) -> case compare exp_re 0 of
             GT -> inf :+ 0
             LT -> 0 :+ 0
             EQ -> nan :+ nan
    | otherwise -> exp (log x * y)
  where
    inf = 1/0
    nan = 0/0

So we must be normally looking at the exp (log x * y) part where exp and log instances of Complex looks like;

exp (x:+y)     =  expx * cos y :+ expx * sin y
                  where expx = exp x
log z          =  log (magnitude z) :+ phase z

then i moved on to the magnitude which is defined like;

magnitude :: (RealFloat a) => Complex a -> a
magnitude (x:+y) =  scaleFloat k
                     (sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
                    where k  = max (exponent x) (exponent y)
                          mk = - k
                          sqr z = z * z

where i am stuck. I would just simply do like sqrt (realPart z ^ 2 + imagPart z ^ 2).

What am i doing wrong..?

Choosing what root to take when dealing with complex number is not a trivial matter. You probably get the principal root. Check Wikipedia to understand the matter better. They have an example with (-8)^(1/3) actually. https://en.m.wikipedia.org/wiki/Nth_root — LudvigH, Sep 18 '17 at 20:58
...but the quick solution to your problem is to write a wrapper around the normal library for positive numbers, and checking that you are taking an odd root of the negative number — LudvigH, Sep 18 '17 at 21:05

K. A. Buhr · Accepted Answer · 2017-09-19T01:34:21.947

Well, the answer you got is a cube root of -8, as you can see by cubing it:

> let x = ((-8):+0)**(1/3)
> x
1.0 :+ 1.732050807568877
> x^3
(-7.9999999999999964) :+ 2.220446049250313e-15
> x*x*x
(-7.9999999999999964) :+ 2.220446049250313e-15
>

It just wasn't the cube root you were looking for.

In fact, there are three cube roots of -8. You can calculate them by calculating the three cube roots of one:

> let j = (-1/2):+sqrt(3)/2   -- a cube root of 1
> let units = [1,j,j*j]       -- that can generate them all

and multiplying them by the value above:

> map (*x) units  -- the cube roots of -8
[1.0 :+ 1.732050807568877,
 (-1.9999999999999996) :+ 1.1102230246251565e-16,
 0.9999999999999997 :+ (-1.7320508075688767)]
> map (^3) $ map (*x) units  -- prove they cube to -8
[(-7.9999999999999964) :+ 2.220446049250313e-15,
 (-7.999999999999995) :+ 1.3322676295501873e-15,
 (-7.999999999999992) :+ 8.881784197001252e-16]
>

As you've determined, the cube root produced by ((-8):+0)**(1/3) is the particular cube root of -8 calculated by:

> exp(log((-8):+0)*(1/3))
1.0 :+ 1.732050807568877
>

It may be worth noting that (1) the definitions of exp and log in Data.Complex, while they may look strange, are the standard mathematical definitions of these functions for complex numbers; and (2) the definition of x ** y as exp(log(x)*y) for complex numbers makes perfect sense mathematically and ensures that x ** y has all the properties you should expect, including the property that cubing ((-8):+0)**(1/3) should give you -8.

As for the definition of magnitude (used in the definition of log), your simpler definition would work too:

> magnitude (-8)
8.0
> sqrt (realPart (-8) ^ 2 + imagPart (-8) ^ 2)
8.0
>

However, the definition of magnitude in Data.Complex has been chosen to provide correct answers in cases where the numbers involved are very large:

> magnitude 1e300
1.0e300
> sqrt (realPart 1e300 ^ 2 + imagPart 1e300 ^ 2)
Infinity
>

Data.Complex package gives strange result for cube root of -8 + 0i ((-8):+0)

1 Answers1