I am trying to divide a double precision floating point number (Double) by 2.0 one hundred times successively in Haskell.
When I get the value 5960464477539.06250, dividing this value by 2.00000 results in 2.9802322387695313e12 (or 2980232238769.53130). Mathematically, the correct result would be 2980232238769.53125, but there is some loss of precision here.
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)
Dividing the values directly in ghci:
Prelude> 5960464477539.06250/2
2.9802322387695313e12
I already read about Data.Scientific, but I don't have access to this module. Without it, is there any way to improve accuracy in this case?
Actually, your number has been computed with full precision; it's just that when printing Doubles, the standard thing to do is to only print as many digits as needed to uniquely identify a Double, and in this case that's fewer than all of them after dividing by 2. You can write your own routine for turning a Double into an exact decimal representation with all digits if you like; you use the usual division+modular reduction trick that's used for showing digits of a normal integer, but with a bit of extra logic for inserting a decimal point when the ratio drops below 1.
But... I've always meant to put together a routine for turning an arbitrary-precision Rational into a sequence of digits in an arbitrary base + information on where the sequence loops, so I took this as an excuse to finally do it. This can also be used to double-check the fact I stated above about full precision being available in your starting code. Here's the data type we'll produce as our end result:
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)
The intention is that d represents the sequence of digits intro d ++ cycle (loop d) (plus some handy auxiliary information). Here's how you can compute it:
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
There are various ways you might want to pretty-print this into a String. Here's one:
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)
We can try it out in ghci:
> pp shows (digits 10 5960464477539.0625) ""
"5.9604644775390625e12"
> pp shows (digits 10 (5960464477539.0625/2)) ""
"2.98023223876953125e12"
Full precision is retained. We can also now verify the first paragraph of the answer by doing all the calculations in Double and then converting to Rational just before printing:
> pp shows (digits 10 (toRational (5960464477539.0625 :: Double))) ""
"5.9604644775390625e12"
> pp shows (digits 10 (toRational (5960464477539.0625/2 :: Double))) ""
"2.98023223876953125e12"
Full precision is retained even for the Double-based calculation.