Seltsames Verhalten von (^) in Haskell

12

Warum gibt GHCi unten eine falsche Antwort?

GHCi

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

Python3

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

UPDATE Ich würde Haskells (^) Funktion wie folgt implementieren.

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

Obwohl meine Version nicht korrekter erscheint als die unten von @WillemVanOnsem bereitgestellte, gibt sie zumindest für diesen speziellen Fall seltsamerweise die richtige Antwort.

Python ist ähnlich.

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))
Zufälliger Typ
quelle
Dies ist ein Fehler der Mantisse. a^24ist ungefähr 2.2437e31, und daher gibt es einen Rundungsfehler, der dies erzeugt.
Willem Van Onsem
Ich verstehe nicht Warum gibt es in GHCi einen Rundungsfehler?
Zufälliger Typ
das hat nichts mit ghci zu tun, es ist einfach so, wie die Gleitkommaeinheit mit Floats umgeht.
Willem Van Onsem
1
Das berechnet, 2.243746917640863e31 - 2.2437469176408626e31was einen kleinen Rundungsfehler hat, der verstärkt wird. Sieht aus wie ein Stornierungsproblem.
Chi
2
Vielleicht verwendet Python einen anderen Algorithmus für die Potenzierung, der in diesem Fall genauer ist? Unabhängig von der verwendeten Sprache weisen Gleitkommaoperationen im Allgemeinen Rundungsfehler auf. Dennoch könnte es interessant sein, die Unterschiede zwischen den beiden Algorithmen zu verstehen.
Chi

Antworten:

14

Kurze Antwort : Es gibt einen Unterschied zwischen (^) :: (Num a, Integral b) => a -> b -> aund (**) :: Floating a => a -> a -> a.

Die (^)Funktion funktioniert nur bei integralen Exponenten. Normalerweise wird ein iterativer Algorithmus verwendet, der jedes Mal prüft, ob die Leistung durch zwei teilbar ist, und die Leistung durch zwei teilt (und wenn nicht teilbar, das Ergebnis mit multipliziert x). Dies bedeutet also, dass für 12insgesamt sechs Multiplikationen durchgeführt werden. Wenn eine Multiplikation einen bestimmten Rundungsfehler aufweist, kann dieser Fehler "explodieren". Wie wir im Quellcode sehen können , ist die (^)Funktion wie folgt implementiert :

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

Die (**)Funktion ist zumindest für Floats und Doubles implementiert, um an der Gleitkommaeinheit zu arbeiten. Wenn wir uns die Implementierung von ansehen (**), sehen wir in der Tat :

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

Dies leitet somit zu der powerFloat# :: Float# -> Float# -> Float#Funktion weiter, die normalerweise vom Compiler mit den entsprechenden FPU-Operationen verknüpft wird.

Wenn wir (**)stattdessen verwenden, erhalten wir auch Null für eine 64-Bit-Gleitkommaeinheit:

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

Wir können zum Beispiel den iterativen Algorithmus in Python implementieren:

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)

Wenn wir dann die gleiche Operation ausführen, erhalte ich lokal:

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

Welches ist der gleiche Wert wie das, was wir (^)in GHCi bekommen.

Willem Van Onsem
quelle
1
Der iterative Algorithmus für (^) in Python gibt diesen Rundungsfehler nicht aus. Ist (*) in Haskell und Python unterschiedlich?
Zufälliger Typ
@Randomdude: Soweit ich weiß, hat die pow(..)Funktion in Python nur einen bestimmten Algorithmus für "int / long", nicht für Floats. Bei Schwimmern wird die Leistung der FPU "zurückgefallen".
Willem Van Onsem
Ich meine, wenn ich die Power-Funktion selbst mit (*) in Python implementiere, genauso wie Haskells Implementierung von (^). Ich benutze keine pow()Funktion.
Zufälliger Typ
2
@Randomdude: Ich habe den Algorithmus in Python implementiert und den gleichen Wert wie in ghc erhalten.
Willem Van Onsem
1
Meine Frage wurde mit meinen Versionen von (^) in Haskell und Python aktualisiert. Gedanken bitte?
Zufälliger Typ