packages feed

splitmix (empty) → 0

raw patch · 8 files changed

+523/−0 lines, 8 filesdep +basedep +base-compatdep +bytestringsetup-changed

Dependencies added: base, base-compat, bytestring, containers, criterion, random, splitmix, tf-random, time

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2017, Oleg Grenrus++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Oleg Grenrus nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,76 @@+# splitmix++Pure Haskell implementation of SplitMix pseudo-random number generator.++## dieharder++> [Dieharder](http://webhome.phy.duke.edu/~rgb/General/dieharder.php) is a random+number generator (rng) testing suite. It is intended to test generators, not+files of possibly random numbers as the latter is a fallacious view of what it+means to be random. Is the number 7 random? If it is generated by a random+process, it might be. If it is made up to serve the purpose of some argument+(like this one) it is not. Perfect random number generators produce "unlikely"+sequences of random numbers – at exactly the right average rate. Testing a rng+is therefore quite subtle.++```+time dieharder-input splitmix | dieharder -a -g 200+```++The test-suite takes around half-an-hour to complete.+All tests are PASSED (occasionally WEAK).++In comparison, built-in [Marsenne Twister](https://en.wikipedia.org/wiki/Mersenne_Twister)+test takes around 15min.++```+time dieharder-input -a+```++## benchmarks++```+benchmarking list/random+time                 96.77 μs   (96.28 μs .. 97.35 μs)+                     1.000 R²   (1.000 R² .. 1.000 R²)+mean                 96.77 μs   (96.35 μs .. 97.63 μs)+std dev              2.001 μs   (1.028 μs .. 3.796 μs)+variance introduced by outliers: 15% (moderately inflated)++benchmarking list/tf-random+time                 60.43 μs   (60.12 μs .. 60.79 μs)+                     1.000 R²   (1.000 R² .. 1.000 R²)+mean                 60.52 μs   (60.26 μs .. 60.91 μs)+std dev              1.120 μs   (780.2 ns .. 1.513 μs)+variance introduced by outliers: 14% (moderately inflated)++benchmarking list/splitmix+time                 16.38 μs   (16.29 μs .. 16.47 μs)+                     1.000 R²   (0.999 R² .. 1.000 R²)+mean                 16.34 μs   (16.25 μs .. 16.51 μs)+std dev              386.8 ns   (201.5 ns .. 669.8 ns)+variance introduced by outliers: 24% (moderately inflated)++benchmarking tree/random+time                 115.4 μs   (108.6 μs .. 126.0 μs)+                     0.942 R²   (0.910 R² .. 0.977 R²)+mean                 116.0 μs   (110.3 μs .. 124.2 μs)+std dev              23.47 μs   (16.60 μs .. 33.95 μs)+variance introduced by outliers: 95% (severely inflated)++benchmarking tree/tf-random+time                 132.4 μs   (132.2 μs .. 132.5 μs)+                     1.000 R²   (1.000 R² .. 1.000 R²)+mean                 132.5 μs   (132.3 μs .. 132.6 μs)+std dev              403.3 ns   (323.0 ns .. 516.5 ns)++benchmarking tree/splitmix+time                 59.80 μs   (59.42 μs .. 60.15 μs)+                     1.000 R²   (0.999 R² .. 1.000 R²)+mean                 59.61 μs   (59.26 μs .. 60.36 μs)+std dev              1.683 μs   (864.5 ns .. 3.206 μs)+variance introduced by outliers: 27% (moderately inflated)+```++Note: the performance can be potentially further improved when GHC gets+[SIMD Support](https://ghc.haskell.org/trac/ghc/wiki/SIMD/Implementation/Status).
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bench/Bench.hs view
@@ -0,0 +1,75 @@+module Main (main) where++import Criterion.Main+import Data.List (unfoldr)+import Data.Word (Word64)++import qualified Data.Tree as T+import qualified System.Random as R+import qualified System.Random.TF as TF+import qualified System.Random.SplitMix as SM++-------------------------------------------------------------------------------+-- List+-------------------------------------------------------------------------------++-- infinite list+genList :: R.RandomGen g => g -> [Int]+genList = unfoldr (Just . R.next)++-- truncated+genListN :: R.RandomGen g => g -> [Int]+genListN = take 2048 . genList++randomList :: Int -> [Int]+randomList = genListN . R.mkStdGen++tfRandomList :: Word64 -> [Int]+tfRandomList w64 = genListN $ TF.seedTFGen (w64, w64, w64, w64)++splitMixList :: Word64 -> [Int]+splitMixList w64 = genListN $ SM.mkSMGen w64++-------------------------------------------------------------------------------+-- Tree+-------------------------------------------------------------------------------++genTree :: R.RandomGen g => g -> T.Tree Int+genTree g = case R.next g of+    ~(i, g') -> T.Node i $ case R.split g' of+        (ga, gb) -> [genTree ga, genTree gb]++genTreeN :: R.RandomGen g => g -> T.Tree Int+genTreeN = cutTree 9 . genTree+  where+    cutTree :: Int -> T.Tree a -> T.Tree a+    cutTree n (T.Node x forest)+        | n <= 0    = T.Node x []+        | otherwise = T.Node x (map (cutTree (n - 1)) forest)++randomTree :: Int -> T.Tree Int+randomTree = genTreeN . R.mkStdGen++tfRandomTree :: Word64 -> T.Tree Int+tfRandomTree w64 = genTreeN $ TF.seedTFGen (w64, w64, w64, w64)++splitMixTree :: Word64 -> T.Tree Int+splitMixTree w64 = genTreeN $ SM.mkSMGen w64++-------------------------------------------------------------------------------+-- Main+-------------------------------------------------------------------------------++main :: IO ()+main = defaultMain+    [ bgroup "list"+        [ bench "random"    $ nf randomList 42+        , bench "tf-random" $ nf tfRandomList 42+        , bench "splitmix"  $ nf splitMixList 42+        ]+    , bgroup "tree"+        [ bench "random"    $ nf randomTree 42+        , bench "tf-random" $ nf tfRandomTree 42+        , bench "splitmix"  $ nf splitMixTree 42+        ]+    ]
+ splitmix.cabal view
@@ -0,0 +1,85 @@+name:		splitmix+version:	0+synopsis:	Fast Splittable PRNG+description:+  Pure Haskell implementation of SplitMix described in+  .+  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 '14). ACM,+  New York, NY, USA, 453-472. DOI:+  <https://doi.org/10.1145/2660193.2660195>+license:	BSD3+license-file:	LICENSE+maintainer:	Oleg Grenrus <oleg.grenrus@iki.fi>+bug-reports:	https://github.com/phadej/splitmix#issues+category:       System+build-type:     Simple+cabal-version:  >= 1.10+tested-with:+  GHC==7.6.3,+  GHC==7.8.4,+  GHC==7.10.3,+  GHC==8.0.2,+  GHC==8.2.1++extra-source-files:+  README.md++library+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   src+  exposed-modules:+    System.Random.SplitMix+  build-depends:+    base   >=4.5     && <4.11,+    time   >=1.4.0.1 && <1.9,+    random >=1.1     && <1.2++  -- dump-core+  -- build-depends: dump-core+  -- ghc-options: -fplugin=DumpCore -fplugin-opt DumpCore:core-html++source-repository head+  type:     git+  location: https://github.com/phadej/splitmix.git++benchmark comparison+  type:             exitcode-stdio-1.0+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   bench+  main-is:          Bench.hs+  build-depends:+    base,+    splitmix,+    random,+    containers >=0.5     && <0.6,+    tf-random  >=0.5     && <0.6,+    criterion  >=1.1.0.0 && <1.3++test-suite montecarlo-pi+  type:             exitcode-stdio-1.0+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   tests+  main-is:          SplitMixPi.hs+  build-depends:+    base,+    splitmix++test-suite dieharder-input+  type:             exitcode-stdio-1.0+  default-language: Haskell2010+  ghc-options:      -Wall+  hs-source-dirs:   tests+  main-is:          Dieharder.hs+  build-depends:+    base,+    splitmix,+    random,+    base-compat >=0.9.3    && <0.10,+    bytestring  >=0.10.4.0 && <0.11,+    tf-random   >=0.5      && <0.6
+ src/System/Random/SplitMix.hs view
@@ -0,0 +1,180 @@+-- |+-- /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>+{-# LANGUAGE Trustworthy #-}+module System.Random.SplitMix (+    SMGen,+    nextWord64,+    nextInt,+    nextDouble,+    splitSMGen,+    -- * Initialisation+    mkSMGen,+    initSMGen,+    newSMGen,+    seedSMGen,+    seedSMGen',+    unseedSMGen,+    ) where++import Data.Bits (popCount, shiftL, shiftR, xor, (.|.))+import Data.IORef (IORef, atomicModifyIORef, newIORef)+import Data.Time.Clock.POSIX (getPOSIXTime)+import Data.Word (Word64, Word32)+import System.CPUTime (getCPUTime, cpuTimePrecision)+import System.IO.Unsafe (unsafePerformIO)++import qualified System.Random as R++-------------------------------------------------------------------------------+-- Generator+-------------------------------------------------------------------------------++-- | SplitMix generator state.+data SMGen = SMGen+    { _seed  :: !Word64+    , _gamma :: !Word64  -- ^ always odd+    }+  deriving Show++-------------------------------------------------------------------------------+-- Operations+-------------------------------------------------------------------------------++-- | Generate a 'Word64'.+nextWord64 :: SMGen -> (Word64, SMGen)+nextWord64 (SMGen seed gamma) = (mix64 seed', SMGen seed' gamma)+  where+    seed' = seed + gamma++-- Here, and in splitSMGen, "inlining" of nextSeed is worth doing.+-- I looked at the Core.+{-+nextWord64 :: SMGen -> (Word64, SMGen)+nextWord64 g0 = (mix64 x, g1)+  where+    (x, g1) = nextSeed g0+-}++-- | Generate an 'Int'.+nextInt :: SMGen -> (Int, SMGen)+nextInt g = case nextWord64 g of+    (w64, g') -> (fromIntegral w64, g')++-- | Generate a 'Double' in @[0, 1)@ range.+nextDouble :: SMGen -> (Double, SMGen)+nextDouble g = case nextWord64 g of+    (w64, g') -> (fromIntegral (w64 `shiftR` 11) * doubleUlp, g')++-- | Split a generator into a two uncorrelated generators.+splitSMGen :: SMGen -> (SMGen, SMGen)+splitSMGen (SMGen seed gamma) =+    (SMGen seed'' gamma, SMGen (mix64 seed') (mixGamma seed''))+  where+    seed'  = seed + gamma+    seed'' = seed' + gamma++{-+splitSMGen :: SMGen -> (SMGen, SMGen)+splitSMGen g0 = (g2, SMGen (mix64 x) (mixGamma y))+  where+    (x, g1) = nextSeed g0+    (y, g2) = nextSeed g1+-}++-------------------------------------------------------------------------------+-- Algorithm+-------------------------------------------------------------------------------++goldenGamma :: Word64+goldenGamma = 0x9e3779b97f4a7c15++doubleUlp :: Double+doubleUlp =  1.0 / fromIntegral (1 `shiftL` 53 :: Word64)++{-+nextSeed :: SMGen -> (Word64, SMGen)+nextSeed (SMGen seed gamma) = (seed', SMGen seed' gamma)+  where+    seed' = seed + gamma+-}++mix64 :: Word64 -> Word64+mix64 z0 =+    let z1 = shiftXorMultiply 33 0xc4ceb9fe1a85ec53 z0+        z2 = shiftXorMultiply 33 0xff51afd7ed558ccd z1+        z3 = shiftXor 33 z2+    in z3++mix64variant13 :: Word64 -> Word64+mix64variant13 z0 =+    let z1 = shiftXorMultiply 30 0xbf58476d1ce4e5b9 z0+        z2 = shiftXorMultiply 27 0x94d049bb133111eb z1+        z3 = shiftXor 31 z2+    in z3++mixGamma :: Word64 -> Word64+mixGamma z0 =+    let z1 = mix64variant13 z0 .|. 1+        n  = popCount (z1 `xor` (z1 `shiftR` 1))+    in if n >= 24+        then z1 `xor` 0xaaaaaaaaaaaaaaaa+        else z1++shiftXor :: Int -> Word64 -> Word64+shiftXor n w = w `xor` (w `shiftR` n)++shiftXorMultiply :: Int -> Word64 -> Word64 -> Word64+shiftXorMultiply n k w = shiftXor n w * k++-------------------------------------------------------------------------------+-- Initialisation+-------------------------------------------------------------------------------++-- | Create 'SMGen' using seed and gamma.+seedSMGen+    :: Word64 -- ^ seed+    -> Word64 -- ^ gamma+    -> SMGen+seedSMGen seed gamma = SMGen seed (gamma .|. 1)++-- | Like 'seedSMGen' but takes a pair.+seedSMGen' :: (Word64, Word64) -> SMGen+seedSMGen' = uncurry seedSMGen++-- | Extract current state of 'SMGen'.+unseedSMGen :: SMGen -> (Word64, Word64)+unseedSMGen (SMGen seed gamma) = (seed, gamma)++-- | Preferred way to deterministically construct 'SMGen'.+mkSMGen :: Word64 -> SMGen+mkSMGen s = SMGen (mix64 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 Word64+mkSeedTime = do+    now <- getPOSIXTime+    cpu <- getCPUTime+    let lo = truncate now :: Word32+        hi = fromIntegral (cpu `div` cpuTimePrecision) :: Word32+    return $ fromIntegral hi `shiftL` 32 .|. fromIntegral lo++-------------------------------------------------------------------------------+-- System.Random+-------------------------------------------------------------------------------++instance R.RandomGen SMGen where+    next = nextInt+    split = splitSMGen
+ tests/Dieharder.hs view
@@ -0,0 +1,47 @@+module Main (main) where++import Prelude+import Prelude.Compat++import Data.List (unfoldr)+import System.Environment (getArgs)+import System.IO (stdout)++import qualified Data.ByteString.Builder as B+import qualified System.Random.SplitMix as SM+import qualified System.Random.TF as TF+import qualified System.Random.TF.Gen as TF+import qualified System.Random.TF.Init as TF++main :: IO ()+main = do+    args <- getArgs+    case args of+        ["splitmix"] -> splitmix+        ["tfrandom"] -> tfrandom+        ["splitmix-tree"] -> splitmixTree+        _ -> return ()++splitmix :: IO ()+splitmix = SM.initSMGen >>= B.hPutBuilder stdout . go+  where+    go :: SM.SMGen -> B.Builder+    go g = case SM.nextWord64 g of +        ~(w64, g') -> B.word64LE w64 `mappend` go g'++splitmixTree :: IO ()+splitmixTree = SM.initSMGen >>= B.hPutBuilder stdout . go+  where+    go :: SM.SMGen -> B.Builder+    go g = case SM.splitSMGen g of+        ~(ga, gb) -> builder 8 ga `mappend` go gb++    builder :: Int -> SM.SMGen -> B.Builder+    builder n = mconcat . take n . map B.word64LE . unfoldr (Just . SM.nextWord64)++tfrandom :: IO ()+tfrandom = TF.initTFGen >>= B.hPutBuilder stdout . go+  where+    go :: TF.TFGen -> B.Builder+    go g = case TF.next g of +        ~(w32, g') -> B.word32LE w32 `mappend` go g'
+ tests/SplitMixPi.hs view
@@ -0,0 +1,28 @@+module Main (main) where++import Data.List (unfoldr, foldl')+import System.Random.SplitMix++doubles :: SMGen -> [Double]+doubles = unfoldr (Just . nextDouble)++monteCarloPi :: SMGen -> Double+monteCarloPi = (4 *) . calc . foldl' accum (P 0 0) . take 10000000 . 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 :: Double)+    print pi'+    print (pi - pi')+