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

实际上，您的数字已经完全精确地计算出来了；只是在打印

Double

时，标准做法是只打印唯一标识

Double

所需的数字，在这种情况下，除以 2 后，数字少于所有数字。你可以写如果你愿意，你自己的例程可以将

Double

转换为所有数字的精确十进制表示；您使用通常的除法+模数缩减技巧，用于显示普通整数的数字，但有一些额外的逻辑用于在比率低于 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

的计算，也能保留完整的精度。

Answer 2

我刚刚以我想要的准确性解决了我的问题（感谢 Daniel Wagner）。

我使用了

Pico

类型 (

Fixed E12

)，分辨率为 10^-12 = .000000000001，来自

Data.Fixed

模块。虽然它不准确，但对于我需要的准确性和易于格式化来说是合理的。这是工作草案：

import Data.Fixed

main :: IO ()
main = do
   n <- readLn :: IO Pico
   reading n 0

roundNew :: Pico -> Integer -> Pico
roundNew n precision = fromIntegral (round (n * 10^precision)) / 10^precision

reading n index = do
   if c == 100
      then return ()
      else do
         putStrLn (show (roundNew n 5))
         reading (n/2.0) (c+1)