diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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). diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -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 + ] + ] diff --git a/splitmix.cabal b/splitmix.cabal new file mode 100644 --- /dev/null +++ b/splitmix.cabal @@ -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: + +license: BSD3 +license-file: LICENSE +maintainer: Oleg Grenrus +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 diff --git a/src/System/Random/SplitMix.hs b/src/System/Random/SplitMix.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/SplitMix.hs @@ -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: +{-# 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 diff --git a/tests/Dieharder.hs b/tests/Dieharder.hs new file mode 100644 --- /dev/null +++ b/tests/Dieharder.hs @@ -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' diff --git a/tests/SplitMixPi.hs b/tests/SplitMixPi.hs new file mode 100644 --- /dev/null +++ b/tests/SplitMixPi.hs @@ -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') +