Haskell中Double值相除时的精度问题

Question

我试图在 Haskell 中将双精度浮点数 ( Double)2.0连续除以一百次。

当我得到值 时5960464477539.06250，将该值除以2.00000结果为2.9802322387695313e12(或2980232238769.53130)。从数学上来说，正确的结果是2980232238769.53125，但这里存在一些精度损失。

import Text.Printf (printf)
import System.Exit (exitSuccess)

main = do
   n <- readLn :: IO Double
   reading n 0

reading n c = do
   if c == 100
      then exitSuccess
      else printf "%.5f\n" n
   reading (n/2.0) (c+1)

将值直接除以ghci：

Prelude> 5960464477539.06250/2
2.9802322387695313e12

我已经阅读过有关 的内容Data.Scientific，但我无权访问该模块。如果没有它，在这种情况下有什么办法可以提高准确性吗？

Answer 1

事实上，你的数字已经被完全精确地计算出来了；只是在打印Doubles 时，标准做法是仅打印唯一标识 a 所需的数字Double，在本例中，该数字少于除以 2 后的所有数字。您可以编写自己的例程来转换 aDouble如果您愿意，可以将其转换为包含所有数字的精确十进制表示形式；您可以使用通常的除法+模数归约技巧，该技巧用于显示普通整数的数字，但需要一些额外的逻辑，以便在比率低于 1 时插入小数点。

但是......我一直想编写一个例程，将任意精度转换Rational为任意基数的数字序列+序列循环位置的信息，所以我以此为借口最终做到了这一点。这也可以用来仔细检查我上面所说的关于起始代码中可用的完整精度的事实。这是我们将生成的最终结果的数据类型：

import Data.Ratio
import Data.Map (Map)
import qualified Data.Map as M

data Sign = Pos | Neg deriving (Bounded, Enum, Eq, Ord, Read, Show)

data Digits = Digits
    { intro :: [Integer]
    , loop :: [Integer]
    , power :: Int
    , sign :: Sign
    , base :: Integer
    } deriving (Eq, Ord, Read, Show)

目的是d表示数字序列intro d ++ cycle (loop d)（加上一些方便的辅助信息）。计算方法如下：

normalizeMagnitude :: Integer -> Rational -> (Rational, Int)
normalizeMagnitude base_ = go where
    base = fromInteger base_
    go x
        | x >= base = succ <$> go (x/base)
        | x < 1 = pred <$> go (x*base)
        | otherwise = (x, 0)

buildLoop :: Integer -> Rational -> (Map Rational (Integer, Rational), Rational)
buildLoop base = go M.empty where
    go seen x = if x `M.member` seen then (seen, x) else go (M.insert x (d, x') seen) x' where
        (d, m) = numerator x `divMod` denominator x
        x' = (m * base) % denominator x

reifyUntil :: Ord a => Map a (b, a) -> a -> a -> [b]
reifyUntil m a0 = go where
    go a = if a0 == a then [] else case M.lookup a m of
        Nothing -> error "never reached sentinel"
        Just (b, a') -> b : go a'

reifyLoop :: Ord a => Map a (b, a) -> a -> [b]
reifyLoop m a0 = case M.lookup a0 m of
    Nothing -> error "not actually a loop"
    Just (b, a) -> b : reifyUntil m a0 a

digits :: Integer -> Rational -> Digits
digits b x = Digits
    { intro = reifyUntil digitMap loopSentinel xNorm
    , loop = reifyLoop digitMap loopSentinel
    , power = p
    , sign = s
    , base = b
    } where
    (xPos, s) = if x < 0 then (-x, Neg) else (x, Pos)
    (xNorm, p) = normalizeMagnitude b xPos
    (digitMap, loopSentinel) = buildLoop b xNorm

您可能想通过多种方式将其漂亮地打印到String. 这是一个：

pp :: (Integer -> ShowS) -> Digits -> ShowS
pp showDigit d = foldr (.) id $ tail [undefined
    , case sign d of
        Pos -> id
        Neg -> ('-':)
    , showDigit x
    , case (xs, xs') of
        ([], [0]) -> id
        _ -> ('.':) . showDigits xs
    , case xs' of
        [0] -> id
        _ -> ('(':) . showDigits xs' . ("...)"++)
    , case (power d, base d) of
        (0, _) -> id
        (_, 10) -> ('e':)
        _ -> ('*':) . shows (base d) . ('^':)
    , case power d of
        0 -> id
        _ -> shows (power d)
    ]
    where
    showDigits = foldr (\digit f -> showDigit digit . f) id
    (x:xs, xs') = case intro d of
        [] -> case loop d of
            [] -> error "invalid digits"
            x:xs -> ([x], xs ++ [x])
        _ -> (intro d, loop d)

我们可以在ghci中尝试一下：

> pp shows (digits 10 5960464477539.0625) ""
"5.9604644775390625e12"
> pp shows (digits 10 (5960464477539.0625/2)) ""
"2.98023223876953125e12"

保留完整精度。Double我们现在还可以通过在打印之前进行所有计算并转换为来验证答案的第一段Rational：

> pp shows (digits 10 (toRational (5960464477539.0625 :: Double))) ""
"5.9604644775390625e12"
> pp shows (digits 10 (toRational (5960464477539.0625/2 :: Double))) ""
"2.98023223876953125e12"

即使基于 - 的计算也能保留完整的精度Double。