Странное поведение (^) в Haskell - PullRequest
Новая
       88

Странное поведение (^) в Haskell

12 голосов
/ 09 января 2020

Почему GHCi дает неправильный ответ ниже?

GHCi 

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

Python3

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

ОБНОВЛЕНИЕ Я бы реализовал функцию Haskell (^) следующим образом. 

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

Хотя моя версия не выглядит более правильной, чем та, что приведена ниже @WillemVanOnsem, она странным образом дает по крайней мере правильный ответ для этого конкретного случая.

Python аналогично. 

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))

1 Ответ

14 голосов
/ 09 января 2020

Краткий ответ : есть разница между (^) :: (Num a, Integral b) => a -> b -> a и (**) :: Floating a => a -> a -> a.

Функция (^) работает только с целыми показателями. Обычно он использует итеративный алгоритм, который каждый раз проверяет, делится ли мощность на два, и делит мощность на два (и, если не делится, умножает результат на x). Таким образом, это означает, что для 12 он выполнит всего шесть умножений. Если умножение имеет определенную ошибку округления, эта ошибка может «взорваться». Как мы видим из исходного кода , функция (^) реализована как : 

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

Функция (**) по крайней мере для Float s и Double s, реализованных для работы с модулем с плавающей запятой. Действительно, если мы посмотрим на реализацию (**), мы увидим: 

instance Floating Float where
    -- &hellip;
    <b>(**) x y = powerFloat x y</b>
    -- &hellip;

Таким образом, это перенаправит на powerFloat# :: Float# -> Float# -> Float# функция, которая обычно будет связана компилятором с соответствующими операциями FPU.

Если вместо этого мы используем (**), мы получим ноль и для 64-битной единицы с плавающей запятой : 

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

Мы можем, например, реализовать итерационный алгоритм в 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)

Если мы затем выполним ту же самую операцию, я получу локально: 

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

Это то же значение, которое мы получаем для (^) в GHCi.

...