pure Splitmix implementation v2
authorCarter Tazio Schonwald <carter.schonwald@gmail.com>
Sun, 15 Jan 2017 15:30:12 +0000 (10:30 -0500)
committerCarter Tazio Schonwald <carter.schonwald@gmail.com>
Sun, 15 Jan 2017 19:43:33 +0000 (14:43 -0500)
random.cabal
src/System/Random/SplitMix/Internal.hs [new file with mode: 0644]

index ac22c68..6f023fd 100644 (file)
@@ -1,10 +1,10 @@
--- Initial random.cabal generated by cabal init.  For further 
+-- Initial random.cabal generated by cabal init.  For further
 -- documentation, see http://haskell.org/cabal/users-guide/
 
 -- The name of the package.
 name:                random
 
--- The package version.  See the Haskell package versioning policy (PVP) 
+-- The package version.  See the Haskell package versioning policy (PVP)
 -- for standards guiding when and how versions should be incremented.
 -- https://wiki.haskell.org/Package_versioning_policy
 -- PVP summary:      +-+------- breaking API changes
@@ -16,7 +16,7 @@ version:             2.0.0.0
 synopsis:            Random number generation library for haskell
 
 -- A longer description of the package.
--- description:         
+-- description:
 
 -- URL for the project homepage or repository.
 homepage:            http://github.com/cartazio/random
@@ -30,18 +30,18 @@ license-file:        LICENSE
 -- The package author(s).
 author:              Carter Tazio Schonwald
 
--- An email address to which users can send suggestions, bug reports, and 
+-- An email address to which users can send suggestions, bug reports, and
 -- patches.
 maintainer:          carter at wellposed dot com
 
 -- A copyright notice.
--- copyright:           
+-- copyright:
 
 category:            Math
 
 build-type:          Simple
 
--- Extra files to be distributed with the package, such as examples or a 
+-- Extra files to be distributed with the package, such as examples or a
 -- README.
 extra-source-files:  ChangeLog.md
 
@@ -51,20 +51,22 @@ cabal-version:       >=1.10
 
 library
   -- Modules exported by the library.
-  -- exposed-modules:     
-  
+  exposed-modules:  System.Random.SplitMix.Internal
+
   -- Modules included in this library but not exported.
-  -- other-modules:       
-  
+  -- other-modules:
+
   -- LANGUAGE extensions used by modules in this package.
-  -- other-extensions:    
-  
+  -- other-extensions:
+
   -- Other library packages from which modules are imported.
-  build-depends:       base >=4.9 && <4.10
-  
+  build-depends:       base >=4.8 && <4.10
+
   -- Directories containing source files.
   hs-source-dirs:      src
-  
+
+  ghc-options: -O2
+
   -- Base language which the package is written in.
   default-language:    Haskell2010
-  
+
diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs
new file mode 100644 (file)
index 0000000..e585586
--- /dev/null
@@ -0,0 +1,156 @@
+{-# LANGUAGE ScopedTypeVariables, BangPatterns, UnboxedTuples #-}
+
+module System.Random.SplitMix.Internal(
+  --mix32,
+  xorShiftR
+  ) where
+
+import qualified  Data.Bits  as DB
+import Data.Bits (xor,(.|.))
+import Data.Word(Word64(..))
+import Data.Functor.Identity
+
+{-# SPECIALIZE popCount :: Word64 -> Word64 #-}
+{-# SPECIALIZE popCount :: Int -> Word64 #-}
+{-# SPECIALIZE popCount :: Word -> Word64 #-}
+popCount :: DB.FiniteBits b => b -> Word64
+popCount  = \ w ->  fromIntegral $ DB.popCount w
+
+
+{-# SPECIALIZE xorShiftR :: Int -> Word64 -> Word64 #-}
+xorShiftR :: DB.FiniteBits  b => Int -> b ->  b
+xorShiftR = \ shift val  ->  val `xor` ( val `DB.unsafeShiftR` shift)
+
+
+xorShiftR33 :: Word64 -> Word64
+xorShiftR33 = \ w -> xorShiftR 33 w
+
+
+firstRoundMix64 :: Word64 -> Word64
+firstRoundMix64 = \ w ->  xorShiftR33 w * 0xff51afd7ed558ccd
+
+secondRoundMix64 :: Word64 -> Word64
+secondRoundMix64 = \ w -> xorShiftR33 w * 0xc4ceb9fe1a85ec53
+
+
+
+mix64variant13 :: Word64 -> Word64
+mix64variant13 = \ w -> xorShiftR 31 $ secondRoundMix64Variant13 $ firstRoundMix64Variant13 w
+
+firstRoundMix64Variant13 :: Word64 -> Word64
+firstRoundMix64Variant13 = \ w -> xorShiftR 30 w * 0xbf58476d1ce4e5b9
+
+
+secondRoundMix64Variant13 :: Word64 -> Word64
+secondRoundMix64Variant13 = \ w -> xorShiftR 27 w * 0x94d049bb133111eb
+
+mix64 :: Word64 -> Word64
+mix64 = \ w -> xorShiftR33 $  secondRoundMix64 $ firstRoundMix64 w
+
+mixGamma :: Word64 -> Word64
+mixGamma = \ w -> runIdentity $!
+  do
+    !mixedGamma <- return $! (mix64variant13 w .|. 1)
+    !bitCount <- return $! popCount $ xorShiftR 1 mixedGamma
+    if bitCount >= 24
+      then return (mixedGamma `xor` 0xaaaaaaaaaaaaaaaa)
+      else return mixedGamma
+
+
+
+data SplitMix64 = SplitMix64 { sm64seed :: {-# UNPACK #-} !Word64
+                              ,sm64Gamma :: {-# UNPACK #-} !Word64 }
+
+
+advanceSplitMix :: SplitMix64 -> SplitMix64
+advanceSplitMix (SplitMix64 sd gamma) = SplitMix64 (sd + gamma) gamma
+
+nextSeedSplitMix :: SplitMix64 -> (# Word64, SplitMix64 #)
+nextSeedSplitMix gen@(SplitMix64 result _) =  newgen `seq` (# result,newgen #)
+  where
+    newgen = advanceSplitMix gen
+
+
+newtype RandomM  a =  RandomM (SplitMix64 -> (# a , SplitMix64 #))
+
+nextWord64SplitMix :: SplitMix64 -> (# Word64 , SplitMix64 #)
+nextWord64SplitMix gen = mixedRes `seq` (# mixedRes , newgen #)
+  where
+    mixedRes = mix64 premixres
+    (#  premixres , newgen  #) = nextSeedSplitMix  gen
+
+splitGeneratorSplitMix :: SplitMix64 -> (# SplitMix64 , SplitMix64 #)
+splitGeneratorSplitMix gen = splitGen `seq`( nextNextGen `seq` (# splitGen , nextNextGen #))
+  where
+    (# splitSeed , nextGen  #) = nextWord64SplitMix gen
+    (# splitPreMixGamma , nextNextGen #) = nextSeedSplitMix nextGen
+    !splitGenGamma = mixGamma splitPreMixGamma
+    !splitGen = SplitMix64 splitSeed splitGenGamma
+
+{-
+
+struct SplitMix64* split_generator(struct SplitMix64* generator) {
+  struct SplitMix64* new_generator = (struct SplitMix64*) malloc(sizeof(struct SplitMix64));
+  new_generator->seed = next_int64(generator);
+  new_generator->gamma = mix_gamma(next_seed(generator));
+  return new_generator;
+}
+
+inline void advance(struct SplitMix64* generator);
+inline uint64_t next_seed(struct SplitMix64* generator);
+
+inline void advance(struct SplitMix64* generator) {
+  generator->seed += generator->gamma;
+}
+
+inline uint64_t next_seed(struct SplitMix64* generator) {
+  uint64_t result = generator->seed;
+  advance(generator);
+  return result;
+}
+
+
+uint64_t next_int64(struct SplitMix64* generator) {
+  return mix64(next_seed(generator));
+}
+
+uint64_t next_bounded_int64(struct SplitMix64* generator, uint64_t bound) {
+  uint64_t threshold = -bound % bound;
+  while (1) {
+    uint64_t r = next_int64(generator);
+    if (r >= threshold) {
+      return r % bound;
+    }
+  }
+}
+
+
+
+struct SplitMix64 {
+  uint64_t seed;
+  uint64_t gamma;
+};
+uint64_t mix_gamma(uint64_t value) {
+  uint64_t mixed_gamma = mix64variant13(value) | 1;
+  int bit_count = pop_count(xor_shift(1, mixed_gamma));
+  if (bit_count >= 24) {
+    return mixed_gamma ^ 0xaaaaaaaaaaaaaaaa;
+  }
+  return mixed_gamma;
+}
+
+uint64_t mix64(uint64_t value) {
+  return xor_shift33(second_round_mix64(first_round_mix64(value)));
+
+inline uint64_t mix64variant13(uint64_t value) {
+  return xor_shift(31, second_round_mix64_variant13(first_round_mix64_variant13(value)));
+
+
+
+inline uint64_t first_round_mix64_variant13(uint64_t value) {
+  return xor_shift(30, value) * 0xbf58476d1ce4e5b9;
+}
+
+inline uint64_t second_round_mix64_variant13(uint64_t value) {
+  return xor_shift(27, value) * 0x94d049bb133111eb;
+-}