haskell浮点计算异常？

Question

使用 ghci 8.6.5

我想计算一个整数输入的平方根，然后将其四舍五入并返回一个整数。

square :: Integer -> Integer
square m = floor $ sqrt $ fromInteger m

它有效。 问题是，对于这个特定的大数字作为输入：

4141414141414141*4141414141414141

我得到了错误的结果。

抛开我的功能，我在 ghci 中测试这个案例：

> sqrt $ fromInteger $ 4141414141414141*4141414141414141
4.1414141414141405e15

错了……对吗？

但很简单

> sqrt $ 4141414141414141*4141414141414141
4.141414141414141e15

这更像是我对计算的期望......

在我的函数中，我必须进行一些类型转换，我认为 fromIntegral 是要走的路。因此，使用它，我的函数为 4141...41 输入给出了错误的结果。

在运行 sqrt 之前，我无法弄清楚 ghci 在类型转换方面的隐含作用。因为 ghci 的转换允许正确计算。

为什么我说这是异常情况：其他号码不会出现问题，例如 5151515151515151 或 3131313131313131 或 4242424242424242 ...

这是一个 Haskell 错误吗？

Answer 1

TLDR

归结为如何将Integer 值转换为无法精确表示的Double。请注意，这不仅可能因为Integer 太大（或太小）而发生，而且Float 和Double 的值通过设计“跳过”整数值，因为它们的大小变得更大。因此，也不是该范围内的每个整数值都可以精确表示。在这种情况下，实现必须根据舍入模式选择一个值。不幸的是，有多个候选人；而你观察到的是 Haskell 选择的候选人会给你一个更差的数字结果。

预期结果

包括 Python 在内的大多数语言都使用所谓的“round-to-nearest-ties-to-even”舍入机制；这是默认的 IEEE754 舍入模式，除非您在兼容处理器中发出浮点相关指令时明确设置舍入模式，否则通常会得到这种模式。在这里使用 Python 作为“参考”，我们得到：

>>> float(long(4141414141414141)*long(4141414141414141))
1.7151311090705027e+31

我还没有尝试过支持所谓的大整数的其他语言，但我希望它们中的大多数都会给你这个结果。

Haskell 如何将Integer 转换为Double

Haskell，然而，使用所谓的截断，或向零舍入。所以你得到：

*Main> (fromIntegral $ 4141414141414141*4141414141414141) :: Double
1.7151311090705025e31

在这种情况下，这是一个“更糟糕”的近似值（参见上面的 Python 生成的值），并且您在原始示例中得到了意想不到的结果。

此时对sqrt 的调用确实是红鲱鱼。

显示代码

一切都源于这段代码：(https://hackage.haskell.org/package/integer-gmp-1.0.2.0/docs/src/GHC.Integer.Type.html#doubleFromInteger)

doubleFromInteger :: Integer -> Double#
doubleFromInteger (S# m#) = int2Double# m#
doubleFromInteger (Jp# bn@(BN# bn#))
    = c_mpn_get_d bn# (sizeofBigNat# bn) 0#
doubleFromInteger (Jn# bn@(BN# bn#))
    = c_mpn_get_d bn# (negateInt# (sizeofBigNat# bn)) 0#

依次调用：(https://github.com/ghc/ghc/blob/master/libraries/integer-gmp/cbits/wrappers.c#L183-L190)：

/* Convert bignum to a `double`, truncating if necessary
 * (i.e. rounding towards zero).
 *
 * sign of mp_size_t argument controls sign of converted double
 */
HsDouble
integer_gmp_mpn_get_d (const mp_limb_t sp[], const mp_size_t sn,
                       const HsInt exponent)
{
...

有目的地表示转换完成向零舍入。

所以，这解释了你得到的行为。

为什么 Haskell 会这样做？

这些都不能解释为什么 Haskell 使用向零舍入进行整数到双精度的转换。我强烈认为它应该使用默认的舍入模式，即round-nearest-ties-to-even。我找不到任何提及这是否是一个有意识的选择，而且它至少与 Python 所做的不同。 （并不是说我认为 Python 是黄金标准，但它确实倾向于让这些事情变得正确。）

我最好的猜测是它只是这样编码的，没有有意识的选择；但也许其他熟悉 Haskell 数字编程历史的人记得更清楚。

做什么

有趣的是，我发现以下讨论可以追溯到 2008 年作为 Python 错误：https://bugs.python.org/issue3166。显然，Python 曾经在这里也做过错误的事情，但他们修复了这种行为。很难追踪确切的历史，但似乎 Haskell 和 Python 都犯了同样的错误； Python 恢复了，但在 Haskell 中没有引起注意。如果这是一个有意识的选择，我想知道为什么。

所以，这就是它的立场。我建议打开一张 GHC 票，这样至少可以正确记录这是“选择”的行为；或者更好，修复它，让它使用默认的舍入模式。

更新：

GHC 票已打开：https://gitlab.haskell.org/ghc/ghc/issues/17231

Answer 2

并非所有Integers 都可以精确地表示为Doubles。对于那些不是，fromInteger 处于需要做出选择的不利位置：它应该返回哪个Double？我在报告中找不到任何讨论在这里做什么的内容，哇！

一个明显的解决方案是返回一个没有小数部分的Double，它表示与存在的任何Double 的原始值的绝对差值最小的整数。不幸的是，这似乎不是 GHC 的fromInteger 做出的决定。

相反，GHC 的选择是返回最大幅度不超过原始数字幅度的Double。所以：

> 17151311090705026844052714160127 :: Double
1.7151311090705025e31
> 17151311090705026844052714160128 :: Double
1.7151311090705027e31

（不要被第二个显示的数字有多短所迷惑：Double 在它上面的行上有整数的确切表示；数字停在那里，因为有足够的唯一标识单个Double。）

为什么这对你很重要？嗯，4141414141414141*4141414141414141 的真实答案是：

> 4141414141414141*4141414141414141
17151311090705026668707274767881

如果fromInteger 将其转换为最接近的Double，如上面的计划（1），它将选择1.7151311090705027e31。但由于它返回的最大 Double 小于上面计划 (2) 中的输入，并且 17151311090705026844052714160128 在技术上更大，因此它返回不太准确的表示 1.7151311090705025e31。

同时，4141414141414141 本身完全可以表示为 Double，因此如果您首先转换为 Double，然后平方，您将得到 Double 选择最接近正确答案的表示的语义，因此计划 (1) 而不是计划 (2)。

这解释了sqrt 输出中的差异：首先在Integer 中进行计算并获得准确答案，然后在最后一秒转换为Double，矛盾的是，它不如转换为Double 准确因为fromInteger 是如何进行转换的，所以立即进行计算并一直进行四舍五入！哎哟。

我怀疑 GHCHQ 会看好修改 fromInteger 以做得更好的补丁；无论如何我知道我会看好它！