Comportamiento extraño de (^) en Haskell

12

¿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))
Tipo al azar
fuente
Este es un error wrt la mantisa. a^24es aproximadamente 2.2437e31, y por lo tanto hay un error de redondeo que produce esto.
Willem Van Onsem
No entiendo. ¿Por qué hay un error de redondeo en GHCi?
Tipo aleatorio el
esto no tiene nada que ver con ghci, es simplemente cómo la unidad de punto flotante maneja los flotadores.
Willem Van Onsem
1
Eso calcula 2.243746917640863e31 - 2.2437469176408626e31que tiene un pequeño error de redondeo que se amplifica. Parece un problema de cancelación.
chi
2
¿Quizás Python usa un algoritmo diferente para la exponenciación, que en este caso es más preciso? En general, sin importar el idioma que utilice, las operaciones de punto flotante exhiben algún error de redondeo. Aún así, podría ser interesante entender las diferencias entre los dos algoritmos.
chi

Respuestas:

14

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.

Willem Van Onsem
fuente
1
El algoritmo iterativo para (^) cuando se implementa en Python no da este error de redondeo. ¿Es (*) diferente en Haskell y Python?
Tipo aleatorio el
@Randomdude: que yo sepa, la pow(..)función en Python solo tiene un cierto algoritmo para "int / long", no para flotantes. Para flotadores, se "retrocederá" en el poder de la FPU.
Willem Van Onsem
Quiero decir, cuando implemento la función de potencia yo mismo usando (*) en Python de la misma manera que la implementación de Haskell de (^). No estoy usando la pow()función.
Tipo aleatorio el
2
@Randomdude: implementé el algoritmo en Python y obtuve el mismo valor que en ghc.
Willem Van Onsem
1
Actualicé mi pregunta con mis versiones de (^) en Haskell y Python. Pensamientos por favor?
Tipo aleatorio el