lundi 20 août 2018

Should Floating-Point RNG be preciser near 0?

The Floating-Point RNG from System.Random looks simple, yet inaccurate to me:

instance Random Double where
  randomR = randomRFloating
  random rng     = 
    case random rng of 
      (x,rng') -> 
          -- We use 53 bits of randomness corresponding to the 53 bit significand:
          ((fromIntegral (mask53 .&. (x::Int64)) :: Double)  
       /  fromIntegral twoto53, rng')
   where 
    twoto53 = (2::Int64) ^ (53::Int64)
    mask53 = twoto53 - 1

Though this RNG indeed produces FP numbers uniformly, there is one thing that I'm doubtful of: There are some numbers in the range that the RNG cannot produce.

Specifically, "too" precise numbers. For example, this RNG can produce (represented as binary IEEE Double-precision FP; sign, exponent, and then mantissa):

0 01111111101 0000000000000000000000000000000000000000000000000000

which is exactly ¼, but cannot produce:

0 01111111101 0000000000000000000000000000000000000000000000000001

because that last 1 (barely) has too much precision.

I doubted if this may happen, so I wrote my own uniform FP RNG:

{-# LANGUAGE ScopedTypeVariables #-}

import Data.Bifunctor
import System.Random

randomFloat1to2 :: (RandomGen g, Random a, RealFloat a) => g -> (a, g) -- Uniformly generates random Float among [1,2)
randomFloat1to2 g = first (1+) (random g)

randomFloatExp :: forall a g. (RandomGen g, Random a, RealFloat a) => Int -> g -> (a, g) -- Uniformly generates random Float among [0, 2^(exp+1))
randomFloatExp exp g = let
    (minexp, _) = floatRange (0 :: a)
    (upperHalf, g') = random g
    in if exp == minexp
        then (0, g') -- Denormal numbers treated as 0
        else if upperHalf
            then first (2^^exp *) (randomFloat1to2 g')
            else randomFloatExp (exp-1) g'

randomFloat :: (RandomGen g, Random a, RealFloat a) => g -> (a, g) -- Uniformly generates random Float among [0,1)
randomFloat = randomFloatExp (-1)

Explanation:

Among the Double numbers in the range [0,1), all numbers in [½,1) have IEEE exponent 01111111110 while the others have lower one. So the RNG makes a coin flip:

If a head comes out, the RNG picks a random number among [½,1) via multiplying ½ and a random number among [1,2). Since the default random effectively picks a random mantissa, we can add 1 to it to make a uniform RNG for the range [1,2).

If not, the RNG makes recursion through [¼,½), [⅛,¼), and so on, until the range is denormal.

Can my version be considered as a better version?




Aucun commentaire:

Enregistrer un commentaire