Haskell中(^)的怪异行为 [英] Weird behavior of (^) in Haskell

问题描述

GHCi为什么在下面给出错误的答案?

Why does GHCi give incorrect answer below?

GHCi

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

Python3

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

更新 我将按以下方式实现Haskell的(^)函数.

UPDATE I would implement Haskell's (^) function as follows.

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在下面提供的版本更正确，但奇怪的是，至少对于这种特殊情况，它给出了正确的答案.

Although my version doesn't appear any more correct than the one provided below by @WillemVanOnsem, it strangely gives the correct answer for this particular case at least.

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

推荐答案

简短答案: (^)函数仅适用于整数指数.通常，它将使用一种迭代算法，该算法每次都会检查幂是否可被二整除，然后将幂除以二(如果不可整，则将结果乘以x).因此，这意味着对于12，它将执行总共六个乘法.如果乘法运算具有一定的舍入误差，则该误差可能会爆炸".正如我们在源代码中看到的， (^)功能实现为:

The (^) function works only on integral exponents. It will normally make use of an iterative algorithm that will each time check if the power is divisible by two, and divide the power by two (and if non-divisible multiply the result with x). This thus means that for 12, it will perform a total of six multiplications. If a multiplication has a certain rounding-off error, that error can "explode". As we can see in the source code, the (^) function is implemented as:

(^) :: (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和Double实现了(**)函数，以在浮点单元上工作.确实，如果我们看一下(**)的实现，我们会看到:

The (**) function is, at least for Floats and Doubles implemented to work on the floating point unit. Indeed, if we take a look at the implementation of (**), we see:

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

因此，这将重定向到

This thus redirect to the powerFloat# :: Float# -> Float# -> Float# function, which will, normally be linked to the corresponding FPU operation(s) by the compiler.

如果我们改用(**)，则对于64位浮点单元，我们也将获得零:

If we use (**) instead, we obtain zero as well for a 64-bit floating point unit:

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

---

例如，我们可以在Python中实现迭代算法:

---

We can for example implement the iterative algorithm in 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)

如果我们随后执行相同的操作，则会在本地获取:

If we then perform the same operation, I get locally:

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

与GHCi中(^)获得的值相同.

Which is the same value as what we get for (^) in GHCi.

这篇关于Haskell中(^)的怪异行为的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！