Dziwne zachowanie (^) w Haskell

12

Dlaczego GHCi podaje poniżej nieprawidłową odpowiedź?

GHCi

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

Python3

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

AKTUALIZACJA Zaimplementowałbym funkcję Haskella (^) w następujący sposób.

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

Chociaż moja wersja nie wydaje się bardziej poprawna niż ta podana poniżej przez @WillemVanOnsem, dziwnie daje co najmniej poprawną odpowiedź dla tego konkretnego przypadku.

Python jest podobny.

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))
Losowy koleś
źródło
To jest błąd wrt mantysy. a^24jest w przybliżeniu 2.2437e31, a zatem występuje błąd zaokrąglania, który to powoduje.
Willem Van Onsem
Nie rozumiem. Dlaczego w GHCi występuje błąd zaokrąglenia?
Losowy koleś
nie ma to nic wspólnego z ghci, jest to po prostu sposób, w jaki jednostka zmiennoprzecinkowa obsługuje zmiennoprzecinkowe.
Willem Van Onsem
1
To oblicza, 2.243746917640863e31 - 2.2437469176408626e31który ma mały błąd zaokrąglania, który zostaje wzmocniony. Wygląda na problem anulowania.
chi
2
Może python używa innego algorytmu do potęgowania, który w tym przypadku jest bardziej precyzyjny? Zasadniczo, bez względu na używany język, operacje zmiennoprzecinkowe wykazują pewne błędy zaokrąglania. Mimo to interesujące może być zrozumienie różnic między tymi dwoma algorytmami.
chi

Odpowiedzi:

14

Krótka odpowiedź : istnieje różnica między (^) :: (Num a, Integral b) => a -> b -> ai (**) :: Floating a => a -> a -> a.

(^)Funkcja działa tylko na integralnych wykładników. Zwykle korzysta z iteracyjnego algorytmu, który za każdym razem sprawdza, czy moc jest podzielna przez dwa, i dzieli moc przez dwa (a jeśli niepodzielna pomnoży wynik x). Oznacza to zatem, że dla 12wykona w sumie sześć multiplikacji. Jeśli mnożenie ma określony błąd zaokrąglania, błąd ten może „wybuchnąć”. Jak widzimy w kodzie źródłowym , (^)funkcja jest implementowana jako :

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

(**)Funkcja jest, przynajmniej dla Floatsi Doubley wdrożony do pracy na jednostkę zmiennoprzecinkową. Rzeczywiście, jeśli spojrzymy na wdrożenie (**), zobaczymy:

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

W ten sposób przekierowuje się do powerFloat# :: Float# -> Float# -> Float#funkcji, która normalnie będzie połączona z odpowiednimi operacjami FPU przez kompilator.

Jeśli użyjemy (**)zamiast tego, otrzymamy również zero dla 64-bitowej jednostki zmiennoprzecinkowej:

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

Możemy na przykład zaimplementować algorytm iteracyjny w Pythonie:

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)

Jeśli następnie wykonamy tę samą operację, otrzymam lokalnie:

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

Jest to ta sama wartość, co otrzymujemy (^)w GHCi.

Willem Van Onsem
źródło
1
Algorytm iteracyjny dla (^) po zaimplementowaniu w Pythonie nie podaje tego błędu zaokrąglania. Czy (*) różni się w Haskell i Python?
Losowy koleś
@Randomdude: o ile wiem, pow(..)funkcja w Pythonie ma tylko pewien algorytm dla „int / long”, a nie dla float. W przypadku pływaków „spadnie” na moc FPU.
Willem Van Onsem
Mam na myśli, kiedy sam implementuję funkcję mocy za pomocą (*) w Pythonie w taki sam sposób, jak implementacja (^) Haskella. Nie używam pow()funkcji.
Losowy koleś
2
@Randomdude: Zaimplementowałem algorytm w Pythonie i uzyskałem taką samą wartość jak w ghc.
Willem Van Onsem
1
Zaktualizowałem moje pytanie o moje wersje (^) w Haskell i Python. Myśli proszę?
Losowy koleś