Weird behavior of (^) in Haskell

Question

Why does GHCi give incorrect answer below?

GHCi

λ> ((-20.24373193905347)^12)^2 - ((-20.24373193905347)^24)
4.503599627370496e15

Python3

>>> ((-20.24373193905347)**12)**2 - ((-20.24373193905347)**24)
0.0

UPDATE I would implement Haskell's (^) function as follows.

powerXY :: Double -> Int -> Double
powerXY x 0 = 1
powerXY x y
    | y < 0 = powerXY (1/x) (-y)
    | otherwise = 
        let z = powerXY x (y `div` 2)
        in  if odd y then z*z*x else z*z

main = do 
    let x = -20.24373193905347
    print $ powerXY (powerXY x 12) 2 - powerXY x 24 -- 0
    print $ ((x^12)^2) - (x ^ 24) -- 4.503599627370496e15

Although my version doesn't appear any more correct than the one provided below by @WillemVanOnsem, it strangely gives the correct answer for this particular case at least.

Python is similar.

def pw(x, y):
    if y < 0:
        return pw(1/x, -y)
    if y == 0:
        return 1
    z = pw(x, y//2)
    if y % 2 == 1:
        return z*z*x
    else:
        return z*z

# prints 0.0
print(pw(pw(-20.24373193905347, 12), 2) - pw(-20.24373193905347, 24))

Answer 1

Short answer: there is a difference between (^) :: (Num a, Integral b) => a -> b -> a and (**) :: Floating a => a -> a -> a.

The (^) function works only on integral exponents. It will normally make use of an iterative algorithm that will each time check if the power is divisible by two, and divide the power by two (and if non-divisible multiply the result with x). This thus means that for 12, it will perform a total of six multiplications. If a multiplication has a certain rounding-off error, that error can "explode". As we can see in the source code, the (^) function is implemented as:

(^) :: (Num a, Integral b) => a -> b -> a
x0 ^ y0 | y0 < 0    = errorWithoutStackTrace "Negative exponent"
        | y0 == 0   = 1
        | otherwise = f x0 y0
    where -- f : x0 ^ y0 = x ^ y
          f x y | even y    = f (x * x) (y `quot` 2)
                | y == 1    = x
                | otherwise = g (x * x) (y `quot` 2) x         -- See Note [Half of y - 1]
          -- g : x0 ^ y0 = (x ^ y) * z
          g x y z | even y = g (x * x) (y `quot` 2) z
                  | y == 1 = x * z
                  | otherwise = g (x * x) (y `quot` 2) (x * z) -- See Note [Half of y - 1]

The (**) function is, at least for Floats and Doubles implemented to work on the floating point unit. Indeed, if we take a look at the implementation of (**), we see:

instance Floating Float where
    -- …
    (**) x y = powerFloat x y
    -- …

This thus redirect to the powerFloat# :: Float# -> Float# -> Float# function, which will, normally be linked to the corresponding FPU operation(s) by the compiler.

If we use (**) instead, we obtain zero as well for a 64-bit floating point unit:

Prelude> (a**12)**2 - a**24
0.0

---

We can for example implement the iterative algorithm in Python:

def pw(x0, y0):
    if y0 < 0:
        raise Error()
    if y0 == 0:
        return 1
    return f(x0, y0)


def f(x, y):
    if (y % 2 == 0):
        return f(x*x, y//2)
    if y == 1:
        return x
    return g(x*x, y // 2, x)


def g(x, y, z):
    if (y % 2 == 0):
        return g(x*x, y//2, z)
    if y == 1:
        return x*z
    return g(x*x, y//2, x*z)

If we then perform the same operation, I get locally:

>>> pw(pw(-20.24373193905347, 12), 2) - pw(-20.24373193905347, 24)
4503599627370496.0

Which is the same value as what we get for (^) in GHCi.