rejection sampling for uniform intervals of smallish integers (those representable...
authorCarter Tazio Schonwald <carter.schonwald@gmail.com>
Mon, 29 May 2017 21:13:55 +0000 (17:13 -0400)
committerCarter Tazio Schonwald <carter.schonwald@gmail.com>
Mon, 29 May 2017 21:13:55 +0000 (17:13 -0400)
random.cabal
src/Data/Distribution/Integers.hs [new file with mode: 0644]

index 597865c..678e561 100644 (file)
@@ -57,6 +57,7 @@ library
     Data.Distribution.FloatingInterval
     Data.Distribution.Normal
     Data.Distribution.Permutation
+    Data.Distribution.Integers
     Data.Random.Utils
 
   c-sources: cmmbits/floatsAndBits.cmm
@@ -71,6 +72,7 @@ library
                       ,ghc-prim
                       ,entropy == 0.3.*
                       ,numeric-extras == 0.1.*
+                      ,primitive >= 0.6
                   -- entropy will later be folded into random, probably
 
   -- Directories containing source files.
diff --git a/src/Data/Distribution/Integers.hs b/src/Data/Distribution/Integers.hs
new file mode 100644 (file)
index 0000000..1f5f122
--- /dev/null
@@ -0,0 +1,30 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+module Data.Distribution.Integers where
+
+import Data.Word(Word64)
+import Data.Bits
+
+--sampleWordRange :: Monad m => (m Word64) -> (Word64,Word64) -> Either
+
+{- | @'sampleWordRangeSimplified' hi@ samples  the closed finite interval @[0,hi]@ -}
+sampleWordRangeSimplified :: forall m . Monad m => m Word64 -> Word64 -> m Word64
+sampleWordRangeSimplified mwd upper
+        | upper + 1 == 0 = mwd
+        -- full 0 ... 2^64-1 range
+        | popCount (upper + 1) == 1 = do  wd <- mwd ; return (wd `mod` (upper + 1))
+        -- power of two range of the form 0 ... 2^k -1, for  0 < k < 64
+        | otherwise = rejectionSampler
+        -- we need to do rejection sampling now!
+    where
+      rejectionSampler :: m Word64
+      rejectionSampler = do awd <- adjustSampleValue
+                            if awd <= upper then return awd
+                              else rejectionSampler
+      nextPower2Log :: Int
+      nextPower2Log =  (64 - countLeadingZeros upper  )
+      adjustSampleValue :: m Word64
+      adjustSampleValue = if nextPower2Log == 64
+                            then mwd
+                            else do wd <- mwd ; return (wd `mod` (bit nextPower2Log))
+
+