Check how many floats in the unit interval are hit by various Word32 -> Float conversion functions

Enumerate.hs

module Main (main) where

import Data.Bit
import Data.Bits
import Data.Vector.Unboxed as V
import Data.Vector.Unboxed.Mutable as MV
import GHC.Float (castFloatToWord32, castWord32ToFloat)
import Data.Word (Word32, Word64)
import Data.Int (Int32)

main :: IO ()
main = do
  putStrLn ""
  putStrLn "random: word32ToFloatInUnitInterval"
  test word32ToFloatInUnitInterval

  putStrLn ""
  putStrLn "mwc-random: wordToFloat"
  test wordToFloat

  putStrLn ""
  putStrLn "downey baseline: wordToFloatByDivision"
  test wordToFloatByDivision

  putStrLn ""
  putStrLn "downey baseline: wordToFloatByDivisionDouble"
  test wordToFloatByDivisionDouble

  putStrLn ""
  putStrLn "java: nextFloat/wordToFloatByDivision24"
  test wordToFloatByDivision24


test :: (Bounded a, Integral a) => (a -> Float) -> IO ()
test f = do
  hits <- generateUpTo maxBound f
  let count = countBits hits
  let full = countBits unitRange
  let overlap = countBits $ zipBits (.&.) hits unitRange

  zeroHit <- indexM hits (fromIntegral (castFloatToWord32 0.0))
  oneHit <- indexM hits (fromIntegral (castFloatToWord32 1.0))

  putStrLn $ "full  = " <> show full
  putStrLn $ "count = " <> show count
  putStrLn $ "percentage hit = " <> show (100 * (fromIntegral count) / (fromIntegral full)) <> "%"
  putStrLn $ "includes 0 = " <> show zeroHit
  putStrLn $ "includes 1 = " <> show oneHit
  putStrLn $ "outside    = " <> show (count - overlap)

false :: Bit
false = zeroBits

true :: Bit
true = complement false

toBit :: Bool -> Bit
toBit True = true
toBit False = false

unitRange :: V.Vector Bit
unitRange = generate (fromIntegral (maxBound :: Word32)) inUnitRange
  where
    inUnitRange :: Int -> Bit
    inUnitRange i = toBit (f >= 0 && f <= 1)
      where
        f = castWord32ToFloat (fromIntegral i)

generateUpTo :: Integral a => a -> (a -> Float) -> IO (V.Vector Bit)
generateUpTo limit f = do
  let go h i
        | i == limit = do
            let j = castFloatToWord32 (f limit)
            () <- MV.unsafeWrite h (fromIntegral j) true
            return h
        | otherwise = do
            let j = castFloatToWord32 (f i)
            () <- MV.unsafeWrite h (fromIntegral j) true
            go h (i + 1)

  hits <- V.unsafeThaw (V.replicate (fromIntegral (maxBound :: Word32)) false)
  go hits 0 >>= V.unsafeFreeze

-- mwc-random
wordToFloat :: Word32 -> Float
wordToFloat x      = (fromIntegral i * m_inv_32) + 0.5 + m_inv_33
    where m_inv_33 = 1.16415321826934814453125e-10
          m_inv_32 = 2.3283064365386962890625e-10
          i        = fromIntegral x :: Int32
{-# INLINE wordToFloat #-}

-- random v1.2 currently
word32ToFloatInUnitInterval :: Word32 -> Float
word32ToFloatInUnitInterval w32 = between1and2 - 1.0
  where
    between1and2 = castWord32ToFloat $ (w32 `unsafeShiftR` 9) .|. 0x3f800000
{-# INLINE word32ToFloatInUnitInterval #-}

-- Downey's baseline
wordToFloatByDivision :: Word32 -> Float
wordToFloatByDivision w = f / m
  where
    f = fromIntegral w :: Float
    m = fromIntegral (maxBound :: Word32) :: Float
{-# INLINE wordToFloatByDivision #-}

-- Downey's baseline using intermediate Double
wordToFloatByDivisionDouble :: Word32 -> Float
wordToFloatByDivisionDouble w = realToFrac d
  where
    f = fromIntegral w :: Double
    m = fromIntegral (maxBound :: Word32) :: Double
    d = f / m :: Double
{-# INLINE wordToFloatByDivisionDouble #-}

-- Java's nextFloat
-- https://docs.oracle.com/javase/8/docs/api/java/util/Random.html#nextFloat--
wordToFloatByDivision24 :: Word32 -> Float
wordToFloatByDivision24 w = f / (fromIntegral m)
  where
    f = fromIntegral (w `unsafeShiftR` 8) :: Float
    m = 1 `unsafeShiftL` 24 :: Word32  
{-# INLINE wordToFloatByDivision24 #-}