Why does random for Float and Double produce exactly 24 or 53 bits? #58
Description
It seems like we could exploit the fact that those types have more precision closer to 0. The algorithm would look like this:
randomDouble :: RandomGen g => g -> (Double, g)
randomDouble = rr where
b :: Word64
b = bit 53
mask = b - 1
r = 1.0 / fromIntegral b
rr g | g `seq` False = undefined
rr g = case randomR (0, mask) g of
(i, g') | testBit i 52 -> seq x (x, g') where
x = r * fromIntegral i
(0, g') -> go0 r 53 g'
(i, g') -> let
cs = countLeadingZeros i - 11
in case randomR (0, bit cs - 1) g' of
(k, g'') -> seq x (x, g'') where
x = r * fromIntegral (unsafeShiftL i cs .|. k) / fromIntegral ((bit cs) :: Word64)
go0 rc sh g | rc `seq` sh `seq` g `seq` False = undefined
-- Stop before hitting denormals, because those are a pain.
go0 rc sh g | sh >= 1022 = (0.0, g)
go0 rc sh g = case randomR (0, mask) g of
(i, g') | testBit i 52 -> seq x (x, g') where
x = rc * fromIntegral i
(0, g') -> go0 (rc * r) (sh + 53) g'
(i, g')
| sh + cs >= 1022 -> (0.0, g')
| otherwise -> case randomR (0, bit cs - 1) g' of
(k, g'') -> seq x (x, g'') where
x = rc * fromIntegral (unsafeShiftL i cs .|. k) / fromIntegral ((bit cs) :: Word64)
where
cs = countLeadingZeros i - 11And then randomRFloating could be defined to take advantage of the increased precision near 0:
randomRFloating :: (Fractional a, Ord a, Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomRFloating = rrf0 where
rrf0 (l, h) g = case compare l h of
LT -> rrf l h g
EQ -> (l, g)
GT -> rrf h l g
rrf l h g | l `seq` h `seq` g `seq` False = undefined
rrf l h g | l >= 0 = case random g of
(coef, g') -> seq x (x, g') where
x = 2.0 * (0.5 * l + coef * (0.5 * h - 0.5 * l))
rrf l h g | h <= 0 = case random g of
(coef, g') -> seq x (x, g') where
x = 2.0 * (0.5 * h + coef * (0.5 * l - 0.5 * h))
-- Here, l < 0 < h. We randomly choose one side and then generate a random number on that side.
rrf l h g = let
rdiv = 1 - toRational h / toRational l
in seq rdiv $ case randomR (0, denominator rdiv - 1) g of
(i, g') | i < numerator rdiv -> go g' where
-- Don't generate 0 on the lower end.
go gc = case random gc of
(r, gc')
| x == 0 = go gc'
| otherwise = (x, gc')
where
x = r * l
(i, g') -> case random g' of
(r, g'') -> seq x (x, g'') where
x = r * h