Haskell 中 (^) 的奇怪行为

Posted 2023-04-14

【中文标题】Haskell 中 (^) 的奇怪行为

【英文标题】：Weird behavior of (^) in Haskell

【发布时间】：2020-01-09 13:07:25

【问题描述】：

为什么 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))

【问题讨论】：

这是一个尾数错误。 a^24 大约是 2.2437e31，因此会产生舍入误差。

我不明白。为什么 GHCi 会出现舍入误差？

这和ghci无关，只是浮点单元如何处理浮点数。

计算2.243746917640863e31 - 2.2437469176408626e31，它有一个小的舍入误差，会被放大。看起来像是取消问题。

也许python使用不同的求幂算法，在这种情况下哪个更精确？一般来说，无论您使用哪种语言，浮点运算都会出现一些舍入误差。不过，了解这两种算法之间的差异可能会很有趣。

【参考方案1】：

简答：(^) :: (Num a, Integral b) => a -> b -> a 和 (**) :: Floating a => a -> a -> a 之间存在差异。

(^) 函数仅适用于整数指数。它通常会使用迭代算法，每次检查幂是否可被 2 整除，然后将幂除以 2（如果不可整除，则将结果与x 相乘）。因此这意味着对于12，它将执行总共六次乘法。如果乘法具有一定的舍入误差，则该误差可能会“爆炸”。正如我们在source code 中看到的，(^) 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]

(**) 函数至少对Floats 和Doubles 实施以在浮点单元上工作。事实上，如果我们看一下(**) 的实现，我们会看到：

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

因此，这将重定向到 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 中为 (^) 得到的值相同。

【讨论】：

在 Python 中实现 (^) 的迭代算法不会产生此舍入错误。 Haskell 和 Python 中的 (*) 是否不同？

@Randomdude：据我所知，Python 中的pow(..) 函数只对“int/long”有一定的算法，对浮点数没有。对于浮点数，它将“回退” FPU 的功能。

我的意思是当我自己在 Python 中使用 (*) 以与 Haskell 的 (^) 实现相同的方式实现幂函数时。我没有使用pow() 函数。

@Randomdude：我已经在 Python 中实现了算法，并获得了与 ghc 中相同的值。

用我在 Haskell 和 Python 中的 (^) 版本更新了我的问题。请问有什么想法吗？