Comportamiento extraño de (^) en Haskell

¿Por qué GHCi da una respuesta incorrecta a continuación?

GHCi

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

Python3

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

ACTUALIZACIÓN Implementaría la función de Haskell (^) de la siguiente manera.

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

Aunque mi versión no parece más correcta que la proporcionada a continuación por @WillemVanOnsem, extrañamente da la respuesta correcta para este caso en particular al menos.

Python es 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))

Respuesta corta : hay una diferencia entre(^) :: (Num a, Integral b) => a -> b -> a y (**) :: Floating a => a -> a -> a.

La (^)función solo funciona en exponentes integrales. Normalmente utilizará un algoritmo iterativo que cada vez verificará si la potencia es divisible por dos, y dividirá la potencia por dos (y si no es divisible multiplica el resultado por x). Esto significa que para 12, realizará un total de seis multiplicaciones. Si una multiplicación tiene un cierto error de redondeo, ese error puede "explotar". Como podemos ver en el código fuente , la (^)función se implementa como :

(^) :: (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]

La (**)función es, al menos para Floats y Doubles implementado para el trabajo en la unidad de coma flotante. De hecho, si echamos un vistazo a la implementación de (**), vemos:

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

Esto, por lo tanto, redirige a la powerFloat# :: Float# -> Float# -> Float# función, que normalmente estará vinculada a las operaciones de FPU correspondientes por el compilador.

Si usamos en su (**)lugar, obtenemos cero también para una unidad de coma flotante de 64 bits:

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

Podemos, por ejemplo, implementar el algoritmo iterativo en 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)

Si luego realizamos la misma operación, obtengo localmente:

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

Cuál es el mismo valor que lo que obtenemos (^)en GHCi.