packages feed

splitmix 0.0.2 → 0.0.3

raw patch · 12 files changed

+787/−64 lines, 12 filesdep +HUnitdep +clockdep ~basedep ~randomdep ~time

Dependencies added: HUnit, clock

Dependency ranges changed: base, random, time

Files

Changelog.md view
@@ -1,3 +1,11 @@+# 0.0.3++- Add `System.Random.SplitMix32` module+- Add `bitmaskWithRejection32` and `bitmaskWithRejection64` module+- Add `nextWord32`, `nextTwoWord32` and `nextFloat`+- Add `random` flag, dropping dependency on `random`+  (breaks things, e.g. `QuickCheck`, when disabled).+ # 0.0.2  - Support back to GHC-7.0
bench/Bench.hs view
@@ -9,6 +9,7 @@ import qualified System.Random.TF as TF import qualified System.Random.TF.Instances as TF import qualified System.Random.SplitMix as SM+import qualified System.Random.SplitMix32 as SM32  ------------------------------------------------------------------------------- -- List@@ -31,6 +32,9 @@ splitMixList :: Word64 -> [Int] splitMixList w64 = genListN $ SM.mkSMGen w64 +splitMix32List :: Word64 -> [Int]+splitMix32List w64 = genListN $ SM32.mkSMGen $ fromIntegral w64+ ------------------------------------------------------------------------------- -- Tree -------------------------------------------------------------------------------@@ -57,6 +61,9 @@ splitMixTree :: Word64 -> T.Tree Int splitMixTree w64 = genTreeN $ SM.mkSMGen w64 +splitMix32Tree :: Word64 -> T.Tree Int+splitMix32Tree w64 = genTreeN $ SM32.mkSMGen $ fromIntegral w64+ ------------------------------------------------------------------------------- -- List Word64 -------------------------------------------------------------------------------@@ -78,6 +85,9 @@ splitMixList64 :: Word64 -> [Word64] splitMixList64 w64 = genListN64 SM.nextWord64 $ SM.mkSMGen w64 +splitMix32List64 :: Word64 -> [Word64]+splitMix32List64 w64 = genListN64 SM32.nextWord64 $ SM32.mkSMGen $ fromIntegral w64+ ------------------------------------------------------------------------------- -- Tree Word64 -------------------------------------------------------------------------------@@ -105,6 +115,9 @@ splitMixTree64 :: Word64 -> T.Tree Word64 splitMixTree64 w64 = genTreeN64 SM.nextWord64 $ SM.mkSMGen w64 +splitMix32Tree64 :: Word64 -> T.Tree Word64+splitMix32Tree64 w64 = genTreeN64 SM32.nextWord64 $ SM32.mkSMGen $ fromIntegral w64+ ------------------------------------------------------------------------------- -- Main -------------------------------------------------------------------------------@@ -112,23 +125,27 @@ main :: IO () main = defaultMain     [ bgroup "list"-        [ bench "random"    $ nf randomList 42-        , bench "tf-random" $ nf tfRandomList 42-        , bench "splitmix"  $ nf splitMixList 42+        [ bench "random"     $ nf randomList 42+        , bench "tf-random"  $ nf tfRandomList 42+        , bench "splitmix"   $ nf splitMixList 42+        , bench "splitmix32" $ nf splitMix32List 42         ]     , bgroup "tree"-        [ bench "random"    $ nf randomTree 42-        , bench "tf-random" $ nf tfRandomTree 42-        , bench "splitmix"  $ nf splitMixTree 42+        [ bench "random"     $ nf randomTree 42+        , bench "tf-random"  $ nf tfRandomTree 42+        , bench "splitmix"   $ nf splitMixTree 42+        , bench "splitmix32" $ nf splitMix32Tree 42         ]     , bgroup "list 64"-        [ bench "random"    $ nf randomList64 42-        , bench "tf-random" $ nf tfRandomList64 42-        , bench "splitmix"  $ nf splitMixList64 42+        [ bench "random"     $ nf randomList64 42+        , bench "tf-random"  $ nf tfRandomList64 42+        , bench "splitmix"   $ nf splitMixList64 42+        , bench "splitmix32" $ nf splitMix32List64 42         ]     , bgroup "tree 64"-        [ bench "random"    $ nf randomTree64 42-        , bench "tf-random" $ nf tfRandomTree64 42-        , bench "splitmix"  $ nf splitMixTree64 42+        [ bench "random"     $ nf randomTree64 42+        , bench "tf-random"  $ nf tfRandomTree64 42+        , bench "splitmix"   $ nf splitMixTree64 42+        , bench "splitmix32" $ nf splitMix32Tree64 42         ]     ]
+ bench/Range.hs view
@@ -0,0 +1,109 @@+-- http://www.pcg-random.org/posts/bounded-rands.html+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP          #-}+module Main where++import Data.Bits+import Data.Bits.Compat+import Data.List        (unfoldr)+import Data.Word        (Word32, Word64)++import qualified System.Random            as R+import qualified System.Random.SplitMix32 as SM++#if defined(__GHCJS__)+#else+import System.Clock (Clock (Monotonic), getTime, toNanoSecs)+import Text.Printf  (printf)+#endif++main :: IO ()+main = do+    gen <- SM.newSMGen++    bench gen (\g h -> R.randomR (0, pred h) g)+    bench gen classicMod+    bench gen intMult+    bench gen bitmaskWithRejection++bench :: g -> (g -> Word32 -> (Word32, g)) -> IO ()+bench gen next = do+    print $ take 70 $ unfoldr (\g -> Just (next g 10)) gen+    clocked $ do+        let x = sumOf next gen+        print x++sumOf :: (g -> Word32 -> (Word32, g)) -> g -> Word32+sumOf next = go 0 2+  where+    go !acc !n g | n > 0xfffff = acc+                 | otherwise    = let (w, g') = next g n in go (acc + w) (succ n) g'++classicMod :: SM.SMGen -> Word32 -> (Word32, SM.SMGen)+classicMod g h =+    let (w32, g') = SM.nextWord32 g in (w32 `mod` h, g')+++-- @+-- uint32_t bounded_rand(rng_t& rng, uint32_t range) {+--     uint32_t x = rng();+--     uint64_t m = uint64_t(x) * uint64_t(range);+--     return m >> 32;+-- }+-- @+--+intMult :: SM.SMGen -> Word32 -> (Word32, SM.SMGen)+intMult g h =+    (fromIntegral $ (fromIntegral w32 * fromIntegral h :: Word64) `shiftR` 32, g')+  where+    (w32, g') = SM.nextWord32 g++-- @+-- uint32_t bounded_rand(rng_t& rng, uint32_t range) {+--     uint32_t mask = ~uint32_t(0);+--     --range;+--     mask >>= __builtin_clz(range|1);+--     uint32_t x;+--     do {+--         x = rng() & mask;+--     } while (x > range);+--     return x;+-- }+-- @@+bitmaskWithRejection :: SM.SMGen -> Word32 -> (Word32, SM.SMGen)+bitmaskWithRejection g0 range = go g0+  where+    mask = complement zeroBits `shiftR` countLeadingZeros (range .|. 1)+    go g = let (x, g') = SM.nextWord32 g+               x' = x .&. mask+           in if x' >= range+              then go g'+              else (x', g')++-------------------------------------------------------------------------------+-- Poor man benchmarking with GHC and GHCJS+-------------------------------------------------------------------------------++clocked :: IO () -> IO ()+#if defined(__GHCJS__)+clocked action = do+    start+    action+    stop++foreign import javascript unsafe+    "console.time('loop');"+    start :: IO ()++foreign import javascript unsafe+    "console.timeEnd('loop');"+    stop :: IO ()+#else+clocked action =  do+    start <- getTime Monotonic+    action+    end <- getTime Monotonic+    printf "loop: %.03fms\n"+        $ fromIntegral (toNanoSecs (end - start))+        / (1e6 :: Double)+#endif
+ bench/SimpleSum.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE CPP #-}+module Main (main) where++import System.Environment (getArgs)+import Data.List (foldl')+import Data.Word (Word32)++import qualified System.Random as R+import qualified System.Random.SplitMix as SM+import qualified System.Random.SplitMix32 as SM32++newGen :: a -> (a -> g) -> IO g -> IO g+#if 0+newGen _ _ new = new+#else+newGen seed mk _ = return (mk seed)+#endif++main :: IO ()+main = do+    putStrLn "Summing randoms..."+    getArgs >>= \args -> case args of+        "splitmix"   : _ -> newGen 33 SM.mkSMGen   SM.newSMGen   >>= \g -> print $ benchSum g SM.nextTwoWord32+        "splitmix32" : _ -> newGen 33 SM32.mkSMGen SM32.newSMGen >>= \g -> print $ benchSum g SM32.nextTwoWord32+        "random"     : _ -> R.newStdGen   >>= \g -> print $ benchSum g randomNextTwoWord32++        -- after Closure Compiler getArgs return [] always?+        -- _ -> newGen 33 SM.mkSMGen   SM.newSMGen   >>= \g -> print $ benchSum g SM.nextTwoWord32+        _ -> newGen 33 SM32.mkSMGen SM32.newSMGen >>= \g -> print $ benchSum g SM32.nextTwoWord32+++benchSum :: g -> (g -> (Word32, Word32, g)) -> Word32+benchSum g next = foldl' (+) 0 $ take 10000000 $ unfoldr2 next g++-- | Infinite unfoldr with two element generator+unfoldr2 :: (s -> (a, a, s)) -> s -> [a]+unfoldr2 f = go where+    go s = let (x, y, s') = f s in x : y : go s'++randomNextTwoWord32 :: R.StdGen -> (Word32, Word32, R.StdGen)+randomNextTwoWord32 s0 = (x, y, s2) where+    (x, s1) = R.random s0+    (y, s2) = R.random s1
splitmix.cabal view
@@ -1,6 +1,6 @@ cabal-version:      >=1.10 name:               splitmix-version:            0.0.2+version:            0.0.3 synopsis:           Fast Splittable PRNG description:   Pure Haskell implementation of SplitMix described in@@ -30,31 +30,49 @@ license-file:       LICENSE maintainer:         Oleg Grenrus <oleg.grenrus@iki.fi> bug-reports:        https://github.com/phadej/splitmix#issues-category:           System+category:           System, Random build-type:         Simple tested-with:-  GHC ==8.6.4 || ==8.4.4 || ==8.2.2 || ==8.0.2 || ==7.10.3 || ==7.8.4 || ==7.6.3 || ==7.4.2 || ==7.2.2 || ==7.0.4+  GHC ==8.8.1 || ==8.6.5 || ==8.4.4 || ==8.2.2 || ==8.0.2 || ==7.10.3 || ==7.8.4 || ==7.6.3 || ==7.4.2 || ==7.2.2 || ==7.0.4  extra-source-files:   README.md   Changelog.md +flag optimised-mixer+  description: Use JavaScript for mix32+  manual:      True+  default:     False++flag random+  description: Providen RandomGen SMGen instance+  manual:      True+  default:     True+ library   default-language: Haskell2010   ghc-options:      -Wall-  hs-source-dirs:   src-  exposed-modules:  System.Random.SplitMix+  hs-source-dirs:   src src-compat+  other-modules:    Data.Bits.Compat+  exposed-modules:+    System.Random.SplitMix+    System.Random.SplitMix32    -- dump-core   -- build-depends: dump-core   -- ghc-options: -fplugin=DumpCore -fplugin-opt DumpCore:core-html    build-depends:-      base     >=4.3     && <4.13+      base     >=4.3     && <4.14     , deepseq  >=1.3.0.0 && <1.5-    , random   >=1.0     && <1.2-    , time     >=1.2.0.3 && <1.9+    , time     >=1.2.0.3 && <1.10 +  if flag(random)+    build-depends: random >=1.0 && <1.2++  if flag(optimised-mixer)+    cpp-options: -DOPTIMISED_MIX32=1+ source-repository head   type:     git   location: https://github.com/phadej/splitmix.git@@ -73,6 +91,41 @@     , splitmix     , tf-random   >=0.5     && <0.6 +benchmark simple-sum+  type:             exitcode-stdio-1.0+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   bench+  main-is:          SimpleSum.hs+  build-depends:+      base+    , random+    , splitmix++benchmark range+  type:             exitcode-stdio-1.0+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   bench src-compat+  main-is:          Range.hs+  other-modules:    Data.Bits.Compat+  build-depends:+      base+    , clock     >=0.8 && <0.9+    , random+    , splitmix++test-suite examples+  type:             exitcode-stdio-1.0+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   tests+  main-is:          Examples.hs+  build-depends:+      base+    , HUnit     ==1.3.1.2 || >=1.6.0.0 && <1.7+    , splitmix+ test-suite montecarlo-pi   type:             exitcode-stdio-1.0   default-language: Haskell2010@@ -83,9 +136,19 @@       base     , splitmix -test-suite splitmix-dieharder+test-suite montecarlo-pi-32   type:             exitcode-stdio-1.0   default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   tests+  main-is:          SplitMixPi32.hs+  build-depends:+      base+    , splitmix++test-suite splitmix-dieharder+  default-language: Haskell2010+  type:             exitcode-stdio-1.0   ghc-options:      -Wall -threaded -rtsopts   hs-source-dirs:   tests   main-is:          Dieharder.hs
+ src-compat/Data/Bits/Compat.hs view
@@ -0,0 +1,45 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-}+module Data.Bits.Compat (+    popCount,+    zeroBits,+    finiteBitSize,+    countLeadingZeros,+    ) where++import Data.Bits++#if !MIN_VERSION_base(4,7,0)+#define FiniteBits Bits+#endif++#if !MIN_VERSION_base(4,5,0)+popCount :: Bits a => a -> Int+popCount = go 0+ where+   go !c 0 = c+   go c w = go (c+1) (w .&. (w - 1)) -- clear the least significant+{-# INLINE popCount #-}+#endif++#if !MIN_VERSION_base(4,7,0)+zeroBits :: Bits a => a+zeroBits = clearBit (bit 0) 0+{-# INLINE zeroBits #-}++finiteBitSize :: Bits a => a -> Int+finiteBitSize = bitSize+{-# INLINE finiteBitSize #-}+#endif++#if !MIN_VERSION_base(4,8,0)+countLeadingZeros :: FiniteBits b => b -> Int+countLeadingZeros x = (w-1) - go (w-1)+  where+    go i | i < 0       = i -- no bit set+         | testBit x i = i+         | otherwise   = go (i-1)++    w = finiteBitSize x+{-# INLINE countLeadingZeros #-}+#endif
src/System/Random/SplitMix.hs view
@@ -1,8 +1,6 @@ -- | -- /SplitMix/ is a splittable pseudorandom number generator (PRNG) that is quite fast. ----- Guy L. Steele, Jr., Doug Lea, and Christine H. Flood. 2014.  Fast splittable pseudorandom number generators. /In Proceedings of the 2014 ACM International Conference on Object Oriented Programming Systems Languages & Applications/ (OOPSLA '13). ACM, New York, NY, USA, 453-472. DOI: <https://doi.org/10.1145/2660193.2660195>--- -- Guy L. Steele, Jr., Doug Lea, and Christine H. Flood. 2014. --  Fast splittable pseudorandom number generators. In Proceedings --  of the 2014 ACM International Conference on Object Oriented@@ -28,17 +26,22 @@ --  but GHC-7.0 and GHC-7.2 have slow implementation, as there --  are no native 'popCount'. ---{-# LANGUAGE CPP #-}-{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP          #-} #if __GLASGOW_HASKELL__ >= 702-{-# LANGUAGE Trustworthy #-}+{-# LANGUAGE Trustworthy  #-} #endif module System.Random.SplitMix (     SMGen,     nextWord64,+    nextWord32,+    nextTwoWord32,     nextInt,     nextDouble,+    nextFloat,     splitSMGen,+    -- * Generation+    bitmaskWithRejection32,+    bitmaskWithRejection64,     -- * Initialisation     mkSMGen,     initSMGen,@@ -49,24 +52,19 @@     ) where  import Control.DeepSeq       (NFData (..))-import Data.Bits             (shiftL, shiftR, xor, (.|.))+import Data.Bits             (complement, shiftL, shiftR, xor, (.&.), (.|.))+import Data.Bits.Compat      (countLeadingZeros, popCount, zeroBits) import Data.IORef            (IORef, atomicModifyIORef, newIORef) import Data.Time.Clock.POSIX (getPOSIXTime) import Data.Word             (Word32, Word64)-import System.CPUTime        (cpuTimePrecision, getCPUTime) import System.IO.Unsafe      (unsafePerformIO) +#ifdef MIN_VERSION_random import qualified System.Random as R+#endif -#if MIN_VERSION_base(4,5,0)-import Data.Bits             (popCount)-#else-import Data.Bits             ((.&.))-popCount :: Word64 -> Int-popCount = go 0- where-   go !c 0 = c-   go c w = go (c+1) (w .&. (w - 1)) -- clear the least significant+#if !__GHCJS__+import System.CPUTime (cpuTimePrecision, getCPUTime) #endif  -- $setup@@ -98,12 +96,12 @@ -- instance Read SMGen where     readsPrec d r =  readParen (d > 10) (\r0 ->-      [ (SMGen seed gamma, r3)-      | ("SMGen", r1) <- lex r0-      , (seed, r2) <- readsPrec 11 r1-      , (gamma, r3) <- readsPrec 11 r2-      , odd gamma-      ]) r+        [ (SMGen seed gamma, r3)+        | ("SMGen", r1) <- lex r0+        , (seed, r2) <- readsPrec 11 r1+        , (gamma, r3) <- readsPrec 11 r2+        , odd gamma+        ]) r  ------------------------------------------------------------------------------- -- Operations@@ -119,6 +117,20 @@   where     seed' = seed + gamma +-- | Generate 'Word32' by truncating 'nextWord64'.+--+-- @since 0.0.3+nextWord32 :: SMGen -> (Word32, SMGen)+nextWord32 g = (fromIntegral w64, g') where+    (w64, g') = nextWord64 g++-- | Generate two 'Word32'.+--+-- @since 0.0.3+nextTwoWord32 :: SMGen -> (Word32, Word32, SMGen)+nextTwoWord32 g = (fromIntegral $ w64 `shiftR` 32, fromIntegral w64, g') where+    (w64, g') = nextWord64 g+ -- | Generate an 'Int'. nextInt :: SMGen -> (Int, SMGen) nextInt g = case nextWord64 g of@@ -133,6 +145,16 @@ nextDouble g = case nextWord64 g of     (w64, g') -> (fromIntegral (w64 `shiftR` 11) * doubleUlp, g') +-- | Generate a 'Float' in @[0, 1)@ range.+--+-- >>> take 8 $ map (printf "%0.3f") $ unfoldr (Just . nextFloat) (mkSMGen 1337) :: [String]+-- ["0.057","0.089","0.237","0.383","0.680","0.320","0.826","0.007"]+--+-- @since 0.0.3+nextFloat :: SMGen -> (Float, SMGen)+nextFloat g = case nextWord32 g of+    (w32, g') -> (fromIntegral (w32 `shiftR` 8) * floatUlp, g')+ -- | Split a generator into a two uncorrelated generators. splitSMGen :: SMGen -> (SMGen, SMGen) splitSMGen (SMGen seed gamma) =@@ -148,13 +170,14 @@ goldenGamma :: Word64 goldenGamma = 0x9e3779b97f4a7c15 +floatUlp :: Float+floatUlp =  1.0 / fromIntegral (1 `shiftL` 24 :: Word32)+ doubleUlp :: Double doubleUlp =  1.0 / fromIntegral (1 `shiftL` 53 :: Word64)  -- Note: in JDK implementations the mix64 and mix64variant13 -- (which is inlined into mixGamma) are swapped.------ I have no idea if swapping them affects statistical properties. mix64 :: Word64 -> Word64 mix64 z0 =    -- MurmurHash3Mixer@@ -168,6 +191,8 @@ mix64variant13 z0 =    -- Better Bit Mixing - Improving on MurmurHash3's 64-bit Finalizer    -- http://zimbry.blogspot.fi/2011/09/better-bit-mixing-improving-on.html+   --+   -- Stafford's Mix13     let z1 = shiftXorMultiply 30 0xbf58476d1ce4e5b9 z0 -- MurmurHash3 mix constants         z2 = shiftXorMultiply 27 0x94d049bb133111eb z1         z3 = shiftXor 31 z2@@ -190,6 +215,41 @@ shiftXorMultiply n k w = shiftXor n w * k  -------------------------------------------------------------------------------+-- Generation+-------------------------------------------------------------------------------++-- | /Bitmask with rejection/ method of generating subrange of 'Word32'.+--+-- @since 0.0.3+bitmaskWithRejection32 :: Word32 -> SMGen -> (Word32, SMGen)+bitmaskWithRejection32 range = go where+    mask = complement zeroBits `shiftR` countLeadingZeros (range .|. 1)+    go g = let (x, g') = nextWord32 g+               x' = x .&. mask+           in if x' >= range+              then go g'+              else (x', g')++-- | /Bitmask with rejection/ method of generating subrange of 'Word64'.+--+-- @bitmaskWithRejection64 w64@ generates random numbers in closed-open+-- range of @[0, w64)@.+--+-- >>> take 20 $ unfoldr (Just . bitmaskWithRejection64 5) (mkSMGen 1337)+-- [3,1,4,1,2,3,1,1,0,3,4,2,3,0,2,3,3,4,1,0]+--+-- @since 0.0.3+bitmaskWithRejection64 :: Word64 -> SMGen -> (Word64, SMGen)+bitmaskWithRejection64 range = go where+    mask = complement zeroBits `shiftR` countLeadingZeros (range .|. 1)+    go g = let (x, g') = nextWord64 g+               x' = x .&. mask+           in if x' >= range+              then go g'+              else (x', g')+++------------------------------------------------------------------------------- -- Initialisation ------------------------------------------------------------------------------- @@ -235,15 +295,21 @@ mkSeedTime :: IO Word64 mkSeedTime = do     now <- getPOSIXTime-    cpu <- getCPUTime     let lo = truncate now :: Word32-        hi = fromIntegral (cpu `div` cpuTimePrecision) :: Word32+#if __GHCJS__+    let hi = lo+#else+    cpu <- getCPUTime+    let hi = fromIntegral (cpu `div` cpuTimePrecision) :: Word32+#endif     return $ fromIntegral hi `shiftL` 32 .|. fromIntegral lo  ------------------------------------------------------------------------------- -- System.Random ------------------------------------------------------------------------------- +#ifdef MIN_VERSION_random instance R.RandomGen SMGen where     next = nextInt     split = splitSMGen+#endif
+ src/System/Random/SplitMix32.hs view
@@ -0,0 +1,311 @@+-- |+-- /SplitMix/ is a splittable pseudorandom number generator (PRNG) that is quite fast.+--+-- This is 32bit variant (original one is 32 bit).+--+-- You __really don't want to use this one__.+--+--  Note: This module supports all GHCs since GHC-7.0.4,+--  but GHC-7.0 and GHC-7.2 have slow implementation, as there+--  are no native 'popCount'.+--+{-# LANGUAGE CPP          #-}+#if __GLASGOW_HASKELL__ >= 702+{-# LANGUAGE Trustworthy  #-}+#endif+module System.Random.SplitMix32 (+    SMGen,+    nextWord32,+    nextWord64,+    nextTwoWord32,+    nextInt,+    nextDouble,+    nextFloat,+    splitSMGen,+    -- * Generation+    bitmaskWithRejection32,+    bitmaskWithRejection64,+    -- * Initialisation+    mkSMGen,+    initSMGen,+    newSMGen,+    seedSMGen,+    seedSMGen',+    unseedSMGen,+    ) where++import Control.DeepSeq       (NFData (..))+import Data.Bits             (complement, shiftL, shiftR, xor, (.&.), (.|.))+import Data.Bits.Compat+       (countLeadingZeros, finiteBitSize, popCount, zeroBits)+import Data.IORef            (IORef, atomicModifyIORef, newIORef)+import Data.Time.Clock.POSIX (getPOSIXTime)+import Data.Word             (Word32, Word64)+import System.IO.Unsafe      (unsafePerformIO)++#ifdef MIN_VERSION_random+import qualified System.Random as R+#endif++#if !__GHCJS__+import System.CPUTime (cpuTimePrecision, getCPUTime)+#endif++-- $setup+-- >>> import Text.Read (readMaybe)+-- >>> import Data.List (unfoldr)+-- >>> import Text.Printf (printf)++-------------------------------------------------------------------------------+-- Generator+-------------------------------------------------------------------------------++-- | SplitMix generator state.+data SMGen = SMGen {-# UNPACK #-} !Word32 {-# UNPACK #-} !Word32 -- seed and gamma; gamma is odd+  deriving Show++instance NFData SMGen where+    rnf (SMGen _ _) = ()++-- |+--+-- >>> readMaybe "SMGen 1 1" :: Maybe SMGen+-- Just (SMGen 1 1)+--+-- >>> readMaybe "SMGen 1 2" :: Maybe SMGen+-- Nothing+--+-- >>> readMaybe (show (mkSMGen 42)) :: Maybe SMGen+-- Just (SMGen 142593372 1604540297)+--+instance Read SMGen where+    readsPrec d r =  readParen (d > 10) (\r0 ->+        [ (SMGen seed gamma, r3)+        | ("SMGen", r1) <- lex r0+        , (seed, r2) <- readsPrec 11 r1+        , (gamma, r3) <- readsPrec 11 r2+        , odd gamma+        ]) r++-------------------------------------------------------------------------------+-- Operations+-------------------------------------------------------------------------------++-- | Generate a 'Word32'.+--+-- >>> take 3 $ map (printf "%x") $ unfoldr (Just . nextWord32) (mkSMGen 1337) :: [String]+-- ["e0cfe722","a6ced0f0","c3a6d889"]+--+nextWord32 :: SMGen -> (Word32, SMGen)+nextWord32 (SMGen seed gamma) = (mix32 seed', SMGen seed' gamma)+  where+    seed' = seed + gamma++-- | Generate a 'Word64', by generating to 'Word32's.+nextWord64 :: SMGen -> (Word64, SMGen)+nextWord64 s0 = (fromIntegral w0 `shiftL` 32 .|. fromIntegral w1,  s2)+  where+    (w0, s1) = nextWord32 s0+    (w1, s2) = nextWord32 s1++-- | Generate two 'Word32'.+nextTwoWord32 :: SMGen -> (Word32, Word32, SMGen)+nextTwoWord32 s0 = (w0, w1, s2) where+    (w0, s1) = nextWord32 s0+    (w1, s2) = nextWord32 s1++-- | Generate an 'Int'.+nextInt :: SMGen -> (Int, SMGen)+nextInt g | isBigInt  = int64+          | otherwise = int32+  where+    int32 = case nextWord32 g of+        (w, g') -> (fromIntegral w, g')+    int64 = case nextWord64 g of+        (w, g') -> (fromIntegral w, g')++isBigInt :: Bool+isBigInt = finiteBitSize (undefined :: Int) > 32++-- | Generate a 'Double' in @[0, 1)@ range.+--+-- >>> take 8 $ map (printf "%0.3f") $ unfoldr (Just . nextDouble) (mkSMGen 1337) :: [String]+-- ["0.878","0.764","0.063","0.845","0.262","0.490","0.176","0.544"]+--+nextDouble :: SMGen -> (Double, SMGen)+nextDouble g = case nextWord64 g of+    (w64, g') -> (fromIntegral (w64 `shiftR` 11) * doubleUlp, g')++-- | Generate a 'Float' in @[0, 1)@ range.+--+-- >>> take 8 $ map (printf "%0.3f") $ unfoldr (Just . nextFloat) (mkSMGen 1337) :: [String]+-- ["0.878","0.652","0.764","0.631","0.063","0.180","0.845","0.645"]+--+nextFloat :: SMGen -> (Float, SMGen)+nextFloat g = case nextWord32 g of+    (w32, g') -> (fromIntegral (w32 `shiftR` 8) * floatUlp, g')++-- | Split a generator into a two uncorrelated generators.+splitSMGen :: SMGen -> (SMGen, SMGen)+splitSMGen (SMGen seed gamma) =+    (SMGen seed'' gamma, SMGen (mix32 seed') (mixGamma seed''))+  where+    seed'  = seed + gamma+    seed'' = seed' + gamma++-------------------------------------------------------------------------------+-- Algorithm+-------------------------------------------------------------------------------++-- | (1 + sqrt 5) / 2 * (2 ^^ bits)+goldenGamma :: Word32+goldenGamma = 0x9e3779b9++floatUlp :: Float+floatUlp =  1.0 / fromIntegral (1 `shiftL` 24 :: Word32)++doubleUlp :: Double+doubleUlp =  1.0 / fromIntegral (1 `shiftL` 53 :: Word64)++#if defined(__GHCJS__) && defined(OPTIMISED_MIX32)+-- JavaScript Foreign Function Interface+-- https://github.com/ghcjs/ghcjs/blob/master/doc/foreign-function-interface.md++foreign import javascript unsafe+    "var x0 = $1 ^ $1 >>> 16; var x1 = x0 & 0xffff; var x2 = (((x0 >>> 16 & 0xffff) * 0x0000ca6b + x1 * 0x000085eb & 0xffff) << 16) + x1 * 0x0000ca6b; var x3 = x2 ^ x2 >>> 13; var x4 = x3 & 0xffff; var x5 = (((x3 >>> 16 & 0xffff) * 0x0000ae35 + x4 * 0x0000c2b2 & 0xffff) << 16) + x4 * 0x0000ae35; $r = (x5 ^ x5 >>> 16) | 0;"+    mix32 :: Word32 -> Word32++foreign import javascript unsafe+    "var x0 = $1 ^ $1 >>> 16; var x1 = x0 & 0xffff; var x2 = (((x0 >>> 16 & 0xffff) * 0x00006ccb + x1 * 0x000069ad & 0xffff) << 16) + x1 * 0x00006ccb; var x3 = x2 ^ x2 >>> 13; var x4 = x3 & 0xffff; var x5 = (((x3 >>> 16 & 0xffff) * 0x0000b5b3 + x4 * 0x0000cd9a & 0xffff) << 16) + x4 * 0x0000b5b3; $r = (x5 ^ x5 >>> 16) | 0;"+    mix32variant13 :: Word32 -> Word32++#else+mix32 :: Word32 -> Word32+mix32 z0 =+   -- MurmurHash3Mixer 32bit+    let z1 = shiftXorMultiply 16 0x85ebca6b z0+        z2 = shiftXorMultiply 13 0xc2b2ae35 z1+        z3 = shiftXor 16 z2+    in z3++-- used only in mixGamma+mix32variant13 :: Word32 -> Word32+mix32variant13 z0 =+   -- See avalanche "executable"+    let z1 = shiftXorMultiply 16 0x69ad6ccb z0+        z2 = shiftXorMultiply 13 0xcd9ab5b3 z1+        z3 = shiftXor 16 z2+    in z3++shiftXor :: Int -> Word32 -> Word32+shiftXor n w = w `xor` (w `shiftR` n)++shiftXorMultiply :: Int -> Word32 -> Word32 -> Word32+shiftXorMultiply n k w = shiftXor n w * k+#endif++mixGamma :: Word32 -> Word32+mixGamma z0 =+    let z1 = mix32variant13 z0 .|. 1             -- force to be odd+        n  = popCount (z1 `xor` (z1 `shiftR` 1))+    -- see: http://www.pcg-random.org/posts/bugs-in-splitmix.html+    -- let's trust the text of the paper, not the code.+    in if n >= 12+        then z1+        else z1 `xor` 0xaaaaaaaa++-------------------------------------------------------------------------------+-- Generation+-------------------------------------------------------------------------------++-- | /Bitmask with rejection/ method of generating subrange of 'Word32'.+bitmaskWithRejection32 :: Word32 -> SMGen -> (Word32, SMGen)+bitmaskWithRejection32 range = go where+    mask = complement zeroBits `shiftR` countLeadingZeros (range .|. 1)+    go g = let (x, g') = nextWord32 g+               x' = x .&. mask+           in if x' >= range+              then go g'+              else (x', g')++-- | /Bitmask with rejection/ method of generating subrange of 'Word64'.+--+-- @bitmaskWithRejection64 w64@ generates random numbers in closed-open+-- range of @[0, w64)@.+--+-- >>> take 20 $ unfoldr (Just . bitmaskWithRejection64 5) (mkSMGen 1337)+-- [0,2,4,2,1,4,2,4,2,2,3,0,3,2,2,2,3,1,2,2]+--+bitmaskWithRejection64 :: Word64 -> SMGen -> (Word64, SMGen)+bitmaskWithRejection64 range = go where+    mask = complement zeroBits `shiftR` countLeadingZeros (range .|. 1)+    go g = let (x, g') = nextWord64 g+               x' = x .&. mask+           in if x' >= range+              then go g'+              else (x', g')++-------------------------------------------------------------------------------+-- Initialisation+-------------------------------------------------------------------------------++-- | Create 'SMGen' using seed and gamma.+--+-- >>> seedSMGen 2 2+-- SMGen 2 3+--+seedSMGen+    :: Word32 -- ^ seed+    -> Word32 -- ^ gamma+    -> SMGen+seedSMGen seed gamma = SMGen seed (gamma .|. 1)++-- | Like 'seedSMGen' but takes a pair.+seedSMGen' :: (Word32, Word32) -> SMGen+seedSMGen' = uncurry seedSMGen++-- | Extract current state of 'SMGen'.+unseedSMGen :: SMGen -> (Word32, Word32)+unseedSMGen (SMGen seed gamma) = (seed, gamma)++-- | Preferred way to deterministically construct 'SMGen'.+--+-- >>> mkSMGen 42+-- SMGen 142593372 1604540297+--+mkSMGen :: Word32 -> SMGen+mkSMGen s = SMGen (mix32 s) (mixGamma (s + goldenGamma))++-- | Initialize 'SMGen' using system time.+initSMGen :: IO SMGen+initSMGen = fmap mkSMGen mkSeedTime++-- | Derive a new generator instance from the global 'SMGen' using 'splitSMGen'.+newSMGen :: IO SMGen+newSMGen = atomicModifyIORef theSMGen splitSMGen++theSMGen :: IORef SMGen+theSMGen = unsafePerformIO $ initSMGen >>= newIORef+{-# NOINLINE theSMGen #-}++mkSeedTime :: IO Word32+mkSeedTime = do+    now <- getPOSIXTime+    let lo = truncate now :: Word32+#if __GHCJS__+    let hi = lo+#else+    cpu <- getCPUTime+    let hi = fromIntegral (cpu `div` cpuTimePrecision) :: Word32+#endif+    return $ fromIntegral hi `shiftL` 32 .|. fromIntegral lo++-------------------------------------------------------------------------------+-- System.Random+-------------------------------------------------------------------------------++#ifdef MIN_VERSION_random+instance R.RandomGen SMGen where+    next = nextInt+    split = splitSMGen+#endif
tests/Dieharder.hs view
@@ -8,6 +8,7 @@  import Control.Concurrent.QSem import Control.DeepSeq         (force)+import Control.Monad           (when) import Data.Bits               (shiftL, (.|.)) import Data.Char               (isSpace) import Data.List               (isInfixOf, unfoldr)@@ -17,7 +18,7 @@ import Foreign.Ptr             (castPtr) import GHC.IO.Exception        (IOErrorType (..), IOException (..)) import System.Environment      (getArgs)-import System.IO               (Handle, hGetContents)+import System.IO               (Handle, hGetContents, stdout) import Text.Printf             (printf)  import qualified Control.Concurrent.Async     as A@@ -27,6 +28,7 @@ import qualified Data.Vector.Storable.Mutable as MSV import qualified System.Process               as Proc import qualified System.Random.SplitMix       as SM+import qualified System.Random.SplitMix32     as SM32 import qualified System.Random.TF             as TF import qualified System.Random.TF.Gen         as TF import qualified System.Random.TF.Init        as TF@@ -37,18 +39,26 @@     if null args     then return ()     else do-        (cmd, runs, conc, seed, test, _help) <- parseArgsIO args $ (,,,,,)+        (cmd, runs, conc, seed, test, raw, _help) <- parseArgsIO args $ (,,,,,,)             <$> arg             <*> optDef "-n" 1             <*> optDef "-j" 1             <*> opt "-s"             <*> opt "-d"+            <*> flag "-r"             <*> flag "-h" +        let run :: RunType g+            run | raw       = runRaw+                | otherwise = runManaged+         case cmd of               "splitmix"      -> do                   g <- maybe SM.initSMGen (return . SM.mkSMGen) seed                   run test runs conc SM.splitSMGen SM.nextWord64 g+              "splitmix32"      -> do+                  g <- maybe SM32.initSMGen (return . SM32.mkSMGen) (fmap fromIntegral seed)+                  run test runs conc SM32.splitSMGen SM32.nextWord64 g               "tfrandom"      -> do                   g <- TF.initTFGen                   run test runs conc TF.split tfNext64 g@@ -64,19 +74,26 @@ -- Dieharder ------------------------------------------------------------------------------- -run :: Maybe Int+type RunType g =+       Maybe Int     -> Int     -> Int     -> (g -> (g, g))     -> (g -> (Word64, g))     -> g-    -> IO ()-run test runs conc split word gen = do+    -> IO () ++runRaw :: RunType g+runRaw _test _runs _conc split word gen =+    generate word split gen stdout++runManaged :: RunType g+runManaged test runs conc split word gen = do     qsem <- newQSem conc      rs <- A.forConcurrently (take runs $ unfoldr (Just . split) gen) $ \g ->         E.bracket_ (waitQSem qsem) (signalQSem qsem) $-            dieharder test (generate word g)+            dieharder test (generate word split g)      case mconcat rs of         Result p w f -> do@@ -84,7 +101,7 @@             printf "PASSED %4d %6.02f%%\n" p (fromIntegral p / total * 100)             printf "WEAK   %4d %6.02f%%\n" w (fromIntegral w / total * 100)             printf "FAILED %4d %6.02f%%\n" f (fromIntegral f / total * 100)-{-# INLINE run #-}+{-# INLINE runManaged #-}  dieharder :: Maybe Int -> (Handle -> IO ()) -> IO Result dieharder test gen = do@@ -144,26 +161,29 @@ size :: Int size = 512 -generate :: forall g. (g -> (Word64, g)) -> g -> Handle -> IO ()-generate word gen0 h = do+generate+    :: forall g. (g -> (Word64, g))+    -> (g -> (g, g))+    -> g -> Handle -> IO ()+generate word split gen0 h = do     vec <- MSV.new size     go gen0 vec   where     go :: g -> MSV.IOVector Word64 -> IO ()     go gen vec = do-        gen' <- write gen vec 0+        let (g1, g2) = split gen+        write g1 vec 0         MSV.unsafeWith vec $ \ptr -> do             bs <- BS.unsafePackCStringLen (castPtr ptr, size * 8)             BS.hPutStr h bs-        go gen' vec+        go g2 vec -    write :: g -> MSV.IOVector Word64 -> Int -> IO g+    write :: g -> MSV.IOVector Word64 -> Int -> IO ()     write !gen !vec !i = do         let (w64, gen') = word gen         MSV.unsafeWrite vec i w64-        if i < size-        then write gen' vec (i + 1)-        else return gen'+        when (i < size) $+            write gen' vec (i + 1) {-# INLINE generate #-}  -------------------------------------------------------------------------------
+ tests/Examples.hs view
@@ -0,0 +1,15 @@+module Main (main) where++import Test.HUnit ((@?=))++import qualified System.Random.SplitMix32 as SM32++main :: IO ()+main = do+    let g = SM32.mkSMGen 42+    show g @?= "SMGen 142593372 1604540297"+    print g++    let (w32, g') = SM32.nextWord32 g+    w32     @?= 1296549791+    show g' @?= "SMGen 1747133669 1604540297"
tests/SplitMixPi.hs view
@@ -7,7 +7,7 @@ doubles = unfoldr (Just . nextDouble)  monteCarloPi :: SMGen -> Double-monteCarloPi = (4 *) . calc . foldl' accum (P 0 0) . take 10000000 . pairs . doubles+monteCarloPi = (4 *) . calc . foldl' accum (P 0 0) . take 50000000 . pairs . doubles   where     calc (P n m) = fromIntegral n / fromIntegral m @@ -25,4 +25,3 @@     print (pi :: Double)     print pi'     print (pi - pi')-
+ tests/SplitMixPi32.hs view
@@ -0,0 +1,27 @@+module Main (main) where++import Data.List (unfoldr, foldl')+import System.Random.SplitMix32++doubles :: SMGen -> [Float]+doubles = unfoldr (Just . nextFloat)++monteCarloPi :: SMGen -> Float+monteCarloPi = (4 *) . calc . foldl' accum (P 0 0) . take 50000000 . pairs . doubles+  where+    calc (P n m) = fromIntegral n / fromIntegral m++    pairs (x : y : xs) = (x, y) : pairs xs+    pairs _ = []++    accum (P n m) (x, y) | x * x + y * y >= 1 = P n (m + 1)+                         | otherwise          = P (n + 1) (m + 1)++data P = P !Int !Int++main :: IO ()+main = do+    pi' <- fmap monteCarloPi newSMGen+    print (pi :: Float)+    print pi'+    print (pi - pi')