diff --git a/random.cabal b/random.cabal index ac22c68160e7136dc65340c65dac43cbae36bbaf..6f023fdee0ecc33540afbbff9e9d1c1e603de549 100644 --- a/random.cabal +++ b/random.cabal @@ -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 index 0000000000000000000000000000000000000000..e585586a85fac6b0b14d6db29c8ae2e58f096b2d --- /dev/null +++ b/src/System/Random/SplitMix/Internal.hs @@ -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; +-}