Comportamento estranho de (^) em Haskell

12

Por que o GHCi fornece respostas incorretas abaixo?

GHCi

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

Python3

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

ATUALIZAÇÃO Eu implementaria a função de Haskell (^) da seguinte maneira.

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

Embora minha versão não pareça mais correta do que a fornecida abaixo por @WillemVanOnsem, estranhamente fornece a resposta correta para este caso em particular, pelo menos.

Python é semelhante.

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))
Cara aleatória
fonte
Este é um erro da mantissa. a^24é aproximadamente 2.2437e31e, portanto, há um erro de arredondamento que produz isso.
Willem Van Onsem
Eu não entendo Por que há um erro de arredondamento no GHCi?
Random du
isso não tem nada a ver com ghci, é simplesmente como a unidade de ponto flutuante lida com flutuadores.
Willem Van Onsem
11
Isso calcula 2.243746917640863e31 - 2.2437469176408626e31que possui um pequeno erro de arredondamento que é amplificado. Parece um problema de cancelamento.
chi
2
Talvez o python use um algoritmo diferente para exponenciação, que neste caso é mais preciso? Em geral, não importa o idioma usado, as operações de ponto flutuante exibem algum erro de arredondamento. Ainda assim, pode ser interessante entender as diferenças entre os dois algoritmos.
chi

Respostas:

14

Resposta curta : há uma diferença entre (^) :: (Num a, Integral b) => a -> b -> ae (**) :: Floating a => a -> a -> a.

A (^)função funciona apenas em expoentes integrais. Normalmente, ele fará uso de um algoritmo iterativo que sempre verificará se a potência é divisível por dois e dividirá a potência por dois (e se não divisível multiplicar o resultado x). Isso significa que 12, para , ele realizará um total de seis multiplicações. Se uma multiplicação tiver um certo erro de arredondamento, esse erro poderá "explodir". Como podemos ver no código fonte , a (^)função é implementada 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]

A (**)função é, pelo menos para Floats e Doubles implementadas para trabalhar na unidade de ponto flutuante. De fato, se dermos uma olhada na implementação de (**), veremos:

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

Assim, ele é redirecionado para a powerFloat# :: Float# -> Float# -> Float#função, que normalmente será vinculada às operações de FPU correspondentes pelo compilador.

Se usarmos (**), obtemos zero também para uma unidade de ponto flutuante de 64 bits:

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

Por exemplo, podemos implementar o algoritmo iterativo em 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)

Se executarmos a mesma operação, eu recebo localmente:

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

Qual é o mesmo valor que obtemos (^)no GHCi.

Willem Van Onsem
fonte
11
O algoritmo iterativo para (^) quando implementado no Python não fornece esse erro de arredondamento. (*) É diferente em Haskell e Python?
Random cara
@ Randomdude: até onde eu sei, a pow(..)função em Python só tem um certo algoritmo para "int / long" s, não para carros alegóricos. Para carros alegóricos, ele "recuará" na potência da FPU.
Willem Van Onsem
Quero dizer, quando eu implemento a função power usando (*) em Python da mesma maneira que a implementação de Haskell (^). Eu não estou usando a pow()função.
Random du
2
@ Randomdude: implementei o algoritmo em Python e obtive o mesmo valor que em ghc.
Willem Van Onsem
11
Atualizei minha pergunta com minhas versões de (^) no Haskell e Python. Pensamentos, por favor?
Random gajo