Comportement étrange de (^) à Haskell

12

Pourquoi GHCi donne-t-il une réponse incorrecte ci-dessous?

GHCi

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

Python3

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

MISE À JOUR J'implémenterais la fonction (^) de Haskell comme suit.

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

Bien que ma version ne semble pas plus correcte que celle fournie ci-dessous par @WillemVanOnsem, elle donne étrangement la bonne réponse pour ce cas particulier au moins.

Python est similaire.

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))
Mec aléatoire
la source
Il s'agit d'une erreur dans la mantisse. a^24est approximativement 2.2437e31, et donc il y a une erreur d'arrondi qui produit cela.
Willem Van Onsem
Je ne comprends pas. Pourquoi y a-t-il une erreur d'arrondi dans GHCi?
Mec aléatoire
cela n'a rien à voir avec ghci, c'est simplement la façon dont l'unité à virgule flottante gère les flotteurs.
Willem Van Onsem
1
Cela calcule 2.243746917640863e31 - 2.2437469176408626e31qui a une petite erreur d'arrondi qui est amplifiée. Ressemble à un problème d'annulation.
chi
2
Peut-être que python utilise un algorithme différent pour l'exponentiation, qui dans ce cas est plus précis? En général, quelle que soit la langue que vous utilisez, les opérations en virgule flottante présentent une erreur d'arrondi. Pourtant, il pourrait être intéressant de comprendre les différences entre les deux algorithmes.
chi

Réponses:

14

Réponse courte : il y a une différence entre (^) :: (Num a, Integral b) => a -> b -> aet (**) :: Floating a => a -> a -> a.

La (^)fonction ne fonctionne que sur les exposants intégraux. Il utilisera normalement un algorithme itératif qui vérifiera à chaque fois si la puissance est divisible par deux et divisera la puissance par deux (et s'il n'est pas divisible, multipliez le résultat par x). Cela signifie donc que pour 12, il effectuera un total de six multiplications. Si une multiplication a une certaine erreur d'arrondi, cette erreur peut "exploser". Comme nous pouvons le voir dans le code source , la (^)fonction est implémentée comme :

(^) :: (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 (**)fonction est, au moins pour Floats et Doubles implémentée pour fonctionner sur l'unité à virgule flottante. En effet, si nous regardons la mise en œuvre de (**), nous voyons:

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

Ceci redirige ainsi vers la powerFloat# :: Float# -> Float# -> Float#fonction, qui sera normalement liée aux opérations FPU correspondantes par le compilateur.

Si nous utilisons à la (**)place, nous obtenons également zéro pour une unité à virgule flottante 64 bits:

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

On peut par exemple implémenter l'algorithme itératif 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 nous effectuons ensuite la même opération, j'obtiens localement:

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

C'est la même valeur que ce que nous obtenons (^)dans GHCi.

Willem Van Onsem
la source
1
L'algorithme itératif pour (^) lorsqu'il est implémenté en Python ne donne pas cette erreur d'arrondi. (*) Est-il différent en Haskell et en Python?
Mec aléatoire
@Randomdude: pour autant que je sache, la pow(..)fonction en Python n'a qu'un certain algorithme pour les "int / longs", pas pour les flottants. Pour les flotteurs, il "retombera" sur la puissance du FPU.
Willem Van Onsem
Je veux dire quand j'implémente la fonction de puissance moi-même en utilisant (*) en Python de la même manière que l'implémentation de Haskell de (^). Je n'utilise pas de pow()fonction.
Mec aléatoire
2
@Randomdude: J'ai implémenté l'algorithme en Python et obtenu la même valeur qu'en ghc.
Willem Van Onsem
1
Mise à jour de ma question avec mes versions de (^) en Haskell et Python. Des pensées s'il vous plait?
Mec aléatoire