Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use SplitMix #2

Closed
wants to merge 9 commits into from
Closed
Show file tree
Hide file tree
Changes from 7 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
143 changes: 21 additions & 122 deletions System/Random.hs
Expand Up @@ -31,10 +31,7 @@
-- instance of 'Random' allows one to generate random values of type
-- 'Float'.
--
-- This implementation uses the Portable Combined Generator of L'Ecuyer
-- ["System.Random\#LEcuyer"] for 32-bit computers, transliterated by
-- Lennart Augustsson. It has a period of roughly 2.30584e18.
--
-- This implementation uses the SplitMix algorithm [1].
-----------------------------------------------------------------------------

#include "MachDeps.h"
Expand Down Expand Up @@ -77,10 +74,12 @@ module System.Random

import Prelude

import Control.Arrow (first)
import Data.Bits
import Data.Int
import Data.Word
import Foreign.C.Types
import qualified System.Random.SplitMix as SM

#ifdef __NHC__
import CPUTime ( getCPUTime )
Expand Down Expand Up @@ -178,23 +177,10 @@ class RandomGen g where
class SplittableGen g where
#endif
-- |The 'split' operation allows one to obtain two distinct random number
-- generators. This is very useful in functional programs (for example, when
-- passing a random number generator down to recursive calls), but very
-- little work has been done on statistically robust implementations of
-- 'split' (["System.Random\#Burton", "System.Random\#Hellekalek"]
-- are the only examples we know of).
-- generators.
split :: g -> (g, g)

{- |
The 'StdGen' instance of 'RandomGen' has a 'genRange' of at least 30 bits.

The result of repeatedly using 'next' should be at least as statistically
robust as the /Minimal Standard Random Number Generator/ described by
["System.Random\#Park", "System.Random\#Carta"].
Until more is known about implementations of 'split', all we require is
that 'split' deliver generators that are (a) not identical and
(b) independently robust in the sense just given.

The 'Show' and 'Read' instances of 'StdGen' provide a primitive way to save the
state of a random number generator.
It is required that @'read' ('show' g) == g@.
Expand All @@ -208,42 +194,18 @@ instance of 'StdGen' has the following properties:
* It guarantees to consume only a finite portion of the string.

* Different argument strings are likely to result in different results.

-}

data StdGen
= StdGen !Int32 !Int32
type StdGen = SM.SMGen

instance RandomGen StdGen where
next = stdNext
genRange _ = stdRange
instance RandomGen SM.SMGen where
next = SM.nextInt

#ifdef ENABLE_SPLITTABLEGEN
instance SplittableGen StdGen where
instance SplittableGen SM.SMGen where
#endif
split = stdSplit

instance Show StdGen where
showsPrec p (StdGen s1 s2) =
showsPrec p s1 .
showChar ' ' .
showsPrec p s2

instance Read StdGen where
readsPrec _p = \ r ->
case try_read r of
r'@[_] -> r'
_ -> [stdFromString r] -- because it shouldn't ever fail.
where
try_read r = do
(s1, r1) <- readDec (dropWhile isSpace r)
(s2, r2) <- readDec (dropWhile isSpace r1)
return (StdGen s1 s2, r2)

{-
If we cannot unravel the StdGen from a string, create
one based on the string given.
-}
split = SM.splitSMGen

stdFromString :: String -> (StdGen, String)
stdFromString s = (mkStdGen num, rest)
where (cs, rest) = splitAt 6 s
Expand All @@ -256,21 +218,10 @@ generator, by mapping an 'Int' into a generator. Again, distinct arguments
should be likely to produce distinct generators.
-}
mkStdGen :: Int -> StdGen -- why not Integer ?
mkStdGen s = mkStdGen32 $ fromIntegral s
mkStdGen s = SM.mkSMGen $ fromIntegral s

{-
From ["System.Random\#LEcuyer"]: "The integer variables s1 and s2 ... must be
initialized to values in the range [1, 2147483562] and [1, 2147483398]
respectively."
-}
mkStdGen32 :: Int32 -> StdGen
mkStdGen32 sMaybeNegative = StdGen (s1+1) (s2+1)
where
-- We want a non-negative number, but we can't just take the abs
-- of sMaybeNegative as -minBound == minBound.
s = sMaybeNegative .&. maxBound
(q, s1) = s `divMod` 2147483562
s2 = q `mod` 2147483398
mkStdGen32 s = SM.mkSMGen $ fromIntegral s

createStdGen :: Integer -> StdGen
createStdGen s = mkStdGen32 $ fromIntegral s
Expand All @@ -286,10 +237,9 @@ Minimal complete definition: 'randomR' and 'random'.
class Random a where
-- | Takes a range /(lo,hi)/ and a random number generator
-- /g/, and returns a random value uniformly distributed in the closed
-- interval /[lo,hi]/, together with a new generator. It is unspecified
-- what happens if /lo>hi/. For continuous types there is no requirement
-- that the values /lo/ and /hi/ are ever produced, but they may be,
-- depending on the implementation and the interval.
-- interval /[lo,hi]/, together with a new generator. For continuous types
-- there is no requirement that the values /lo/ and /hi/ are ever produced,
-- but they may be, depending on the implementation and the interval.
randomR :: RandomGen g => (a,a) -> g -> (a,g)

-- | The same as 'randomR', but using a default range determined by the type:
Expand Down Expand Up @@ -448,12 +398,6 @@ instance Random CDouble where
-- random rng = case random rng of
-- (x,rng') -> (realToFrac (x::Double), rng')

mkStdRNG :: Integer -> IO StdGen
mkStdRNG o = do
ct <- getCPUTime
(sec, psec) <- getTime
return (createStdGen (sec * 12345 + psec + ct + o))

randomBounded :: (RandomGen g, Random a, Bounded a) => g -> (a, g)
randomBounded = randomR (minBound, maxBound)

Expand Down Expand Up @@ -510,39 +454,6 @@ randomIvalDouble (l,h) fromDouble rng
int32Count :: Integer
int32Count = toInteger (maxBound::Int32) - toInteger (minBound::Int32) + 1 -- GHC ticket #3982

stdRange :: (Int,Int)
stdRange = (1, 2147483562)

stdNext :: StdGen -> (Int, StdGen)
-- Returns values in the range stdRange
stdNext (StdGen s1 s2) = (fromIntegral z', StdGen s1'' s2'')
where z' = if z < 1 then z + 2147483562 else z
z = s1'' - s2''

k = s1 `quot` 53668
s1' = 40014 * (s1 - k * 53668) - k * 12211
s1'' = if s1' < 0 then s1' + 2147483563 else s1'

k' = s2 `quot` 52774
s2' = 40692 * (s2 - k' * 52774) - k' * 3791
s2'' = if s2' < 0 then s2' + 2147483399 else s2'

stdSplit :: StdGen -> (StdGen, StdGen)
stdSplit std@(StdGen s1 s2)
= (left, right)
where
-- no statistical foundation for this!
left = StdGen new_s1 t2
right = StdGen t1 new_s2

new_s1 | s1 == 2147483562 = 1
| otherwise = s1 + 1

new_s2 | s2 == 1 = 2147483398
| otherwise = s2 - 1

StdGen t1 t2 = snd (next std)

-- The global random number generator

{- $globalrng #globalrng#
Expand All @@ -564,7 +475,7 @@ getStdGen = readIORef theStdGen

theStdGen :: IORef StdGen
theStdGen = unsafePerformIO $ do
rng <- mkStdRNG 0
rng <- SM.initSMGen
curiousleo marked this conversation as resolved.
Show resolved Hide resolved
newIORef rng

-- |Applies 'split' to the current global random generator,
Expand All @@ -588,22 +499,10 @@ getStdRandom f = atomicModifyIORef' theStdGen (swap . f)

{- $references

1. FW #Burton# Burton and RL Page, /Distributed random number generation/,
Journal of Functional Programming, 2(2):203-212, April 1992.

2. SK #Park# Park, and KW Miller, /Random number generators -
good ones are hard to find/, Comm ACM 31(10), Oct 1988, pp1192-1201.

3. DG #Carta# Carta, /Two fast implementations of the minimal standard
random number generator/, Comm ACM, 33(1), Jan 1990, pp87-88.

4. P #Hellekalek# Hellekalek, /Don\'t trust parallel Monte Carlo/,
Department of Mathematics, University of Salzburg,
<http://random.mat.sbg.ac.at/~peter/pads98.ps>, 1998.

5. Pierre #LEcuyer# L'Ecuyer, /Efficient and portable combined random
number generators/, Comm ACM, 31(6), Jun 1988, pp742-749.

The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
1. 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

-}
2 changes: 2 additions & 0 deletions cabal.project
@@ -0,0 +1,2 @@
packages: .
constraints: splitmix -random
2 changes: 1 addition & 1 deletion random.cabal
Expand Up @@ -36,7 +36,7 @@ Library
System.Random
extensions: CPP
GHC-Options: -O2
build-depends: base >= 3 && < 5, time
build-depends: base >= 3 && < 5, time, splitmix

source-repository head
type: git
Expand Down
8 changes: 8 additions & 0 deletions shell.nix
@@ -0,0 +1,8 @@
with import <nixpkgs> {};

mkShell {
buildInputs = [
cabal-install
haskell.compiler.ghc822
];
}