在Haskell中对双精度值进行除法运算时出现精度问题

tpxzln5u  于 2023-03-19  发布在  其他
关注(0)|答案(2)|浏览(177)

我尝试在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的文章,但是我没有访问这个模块的权限。如果没有它,在这种情况下有什么方法可以提高准确性吗?

ndasle7k

ndasle7k1#

实际上,你的数字已经被精确地计算出来了;只是在打印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的计算,也会保留完全精度。

0sgqnhkj

0sgqnhkj2#

我只是用我想要的准确性解决了我的问题(感谢丹尼尔·瓦格纳)。
我使用了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)

相关问题