Why does random for Float and Double produce exactly 24 or 53 bits?

[Zemyla]

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 - 11

And 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

[cartazio]

its a bug/problem in current random

[cartazio]

a better algorithm with decent perf can be found here https://github.com/cartazio/old-random/blob/v1.3.0/src/Data/Distribution/FloatingInterval.hs

[cartazio]

you raised a very good point @Zemyla :)

[cartazio]

i've been in a rabbit hole on some ghc patching experiments for vector, i'll dust this stuff off post haste :)

[Bodigrim]

@Zemyla could you please compare your algorithm to random-1.2.0? I believe we generate more than 24/53 bits now.

random/src/System/Random/Internal.hs

Lines 771 to 790 in 11464aa

	-- | See [Floating point number caveats](System-Random-Stateful.html#fpcaveats).
	instance UniformRange Double where
	uniformRM (l, h) g
	| l == h = return l
	| otherwise = do
	x <- uniformDouble01M g
	return $ x * l + (1 -x) * h
	-- | Generates uniformly distributed 'Double' in the range \([0, 1]\).
	-- Numbers are generated by generating uniform 'Word64' and dividing
	-- it by \(2^{64}\). It's used to implement 'UniformR' instance for
	-- 'Double'.
	--
	-- @since 1.2.0
	uniformDouble01M :: StatefulGen g m => g -> m Double
	uniformDouble01M g = do
	w64 <- uniformWord64 g
	return $ fromIntegral w64 / m
	where
	m = fromIntegral (maxBound :: Word64) :: Double

This approach does not cover denormalized floats (@curiousleo correct me, if I'm wrong, please), which is discussed in #53, but AFAIU yours does not as well.

[curiousleo]

I don't claim to fully understand @Zemyla's code, but my impression is that it is based on ideas similar to http://allendowney.com/research/rand/ -- the countLeadingZeros here can be used to read a uniform bitstring as a geometric distribution, which is how the exponent of a uniformly random IEEE float ought to be generated.

I would like to see something along those lines in a future version of random, either as the default way to generate floating point numbers or as an optional "more accurate but slower" method, depending on how it performs.

However, code like this needs to be tested thoroughly. In addition, while generating numbers in the unit interval is nice, translating the unit interval into an arbitrary target interval via addition and multiplication leads to a loss in precision. So ideally, the uniform floating point generation method would be able to directly generate numbers in an arbitrary interval.

This is what https://gitlab.com/christoph-conrads/rademacher-fpl/ does. Note how thoroughly tested it is - which is necessary, because the code gets pretty complex.

@Zemyla if this is something you're interested in, you may also be interested in the discussions we had in the run-up to the 1.2.0 release: idontgetoutmuch#113 (comment) (and the whole surrounding issue) as well as idontgetoutmuch#105. I am experimenting with a Haskell implementation of a "uniform floating point number in an arbitrary range generator" here: https://github.com/curiousleo/random-float.