From a967cb08d37da9cd04d6e54202a945bcaae89c05 Mon Sep 17 00:00:00 2001 From: Herbert Valerio Riedel Date: Sun, 21 Dec 2014 10:43:15 +0100 Subject: [PATCH 01/31] M-x delete-trailing-whitespace & M-x untabify --- System/Random.hs | 194 +++++++++++++++++++++++------------------------ 1 file changed, 97 insertions(+), 97 deletions(-) diff --git a/System/Random.hs b/System/Random.hs index ab7727405..dfd20883b 100644 --- a/System/Random.hs +++ b/System/Random.hs @@ -7,7 +7,7 @@ -- Module : System.Random -- Copyright : (c) The University of Glasgow 2001 -- License : BSD-style (see the file LICENSE in the 'random' repository) --- +-- -- Maintainer : libraries@haskell.org -- Stability : stable -- Portability : portable @@ -18,7 +18,7 @@ -- or to get different results on each run by using the system-initialised -- generator or by supplying a seed from some other source. -- --- The library is split into two layers: +-- The library is split into two layers: -- -- * A core /random number generator/ provides a supply of bits. -- The class 'RandomGen' provides a common interface to such generators. @@ -40,40 +40,40 @@ #include "MachDeps.h" module System.Random - ( + ( - -- $intro + -- $intro - -- * Random number generators + -- * Random number generators #ifdef ENABLE_SPLITTABLEGEN - RandomGen(next, genRange) - , SplittableGen(split) + RandomGen(next, genRange) + , SplittableGen(split) #else - RandomGen(next, genRange, split) + RandomGen(next, genRange, split) #endif - -- ** Standard random number generators - , StdGen - , mkStdGen + -- ** Standard random number generators + , StdGen + , mkStdGen - -- ** The global random number generator + -- ** The global random number generator - -- $globalrng + -- $globalrng - , getStdRandom - , getStdGen - , setStdGen - , newStdGen + , getStdRandom + , getStdGen + , setStdGen + , newStdGen - -- * Random values of various types - , Random ( random, randomR, - randoms, randomRs, - randomIO, randomRIO ) + -- * Random values of various types + , Random ( random, randomR, + randoms, randomRs, + randomIO, randomRIO ) - -- * References - -- $references + -- * References + -- $references - ) where + ) where import Prelude @@ -83,15 +83,15 @@ import Data.Word import Foreign.C.Types #ifdef __NHC__ -import CPUTime ( getCPUTime ) +import CPUTime ( getCPUTime ) import Foreign.Ptr ( Ptr, nullPtr ) -import Foreign.C ( CTime, CUInt ) +import Foreign.C ( CTime, CUInt ) #else -import System.CPUTime ( getCPUTime ) -import Data.Time ( getCurrentTime, UTCTime(..) ) +import System.CPUTime ( getCPUTime ) +import Data.Time ( getCurrentTime, UTCTime(..) ) import Data.Ratio ( numerator, denominator ) #endif -import Data.Char ( isSpace, chr, ord ) +import Data.Char ( isSpace, chr, ord ) import System.IO.Unsafe ( unsafePerformIO ) import Data.IORef ( IORef, newIORef, readIORef, writeIORef ) #if MIN_VERSION_base (4,6,0) @@ -99,7 +99,7 @@ import Data.IORef ( atomicModifyIORef' ) #else import Data.IORef ( atomicModifyIORef ) #endif -import Numeric ( readDec ) +import Numeric ( readDec ) #ifdef __GLASGOW_HASKELL__ import GHC.Exts ( build ) @@ -201,17 +201,17 @@ It is required that @'read' ('show' g) == g@. In addition, 'reads' may be used to map an arbitrary string (not necessarily one produced by 'show') onto a value of type 'StdGen'. In general, the 'Read' -instance of 'StdGen' has the following properties: +instance of 'StdGen' has the following properties: -* It guarantees to succeed on any string. +* It guarantees to succeed on any string. -* It guarantees to consume only a finite portion of the string. +* It guarantees to consume only a finite portion of the string. * Different argument strings are likely to result in different results. -} -data StdGen +data StdGen = StdGen !Int32 !Int32 instance RandomGen StdGen where @@ -224,8 +224,8 @@ instance SplittableGen StdGen where split = stdSplit instance Show StdGen where - showsPrec p (StdGen s1 s2) = - showsPrec p s1 . + showsPrec p (StdGen s1 s2) = + showsPrec p s1 . showChar ' ' . showsPrec p s2 @@ -234,11 +234,11 @@ instance Read StdGen where case try_read r of r'@[_] -> r' _ -> [stdFromString r] -- because it shouldn't ever fail. - where + where try_read r = do (s1, r1) <- readDec (dropWhile isSpace r) - (s2, r2) <- readDec (dropWhile isSpace r1) - return (StdGen s1 s2, r2) + (s2, r2) <- readDec (dropWhile isSpace r1) + return (StdGen s1 s2, r2) {- If we cannot unravel the StdGen from a string, create @@ -246,7 +246,7 @@ instance Read StdGen where -} stdFromString :: String -> (StdGen, String) stdFromString s = (mkStdGen num, rest) - where (cs, rest) = splitAt 6 s + where (cs, rest) = splitAt 6 s num = foldl (\a x -> x + 3 * a) 1 (map ord cs) @@ -266,11 +266,11 @@ 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 + -- 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 createStdGen :: Integer -> StdGen createStdGen s = mkStdGen32 $ fromIntegral s @@ -323,7 +323,7 @@ class Random a where -- | A variant of 'random' that uses the global random number generator -- (see "System.Random#globalrng"). randomIO :: IO a - randomIO = getStdRandom random + randomIO = getStdRandom random -- | Produce an infinite list-equivalent of random values. {-# INLINE buildRandoms #-} @@ -340,7 +340,7 @@ buildRandoms cons rand = go instance Random Integer where randomR ival g = randomIvalInteger ival g - random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g + random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g instance Random Int where randomR = randomIvalIntegral; random = randomBounded instance Random Int8 where randomR = randomIvalIntegral; random = randomBounded @@ -378,13 +378,13 @@ instance Random CIntMax where randomR = randomIvalIntegral; random = randomBo instance Random CUIntMax where randomR = randomIvalIntegral; random = randomBounded instance Random Char where - randomR (a,b) g = + randomR (a,b) g = case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of (x,g') -> (chr x, g') - random g = randomR (minBound,maxBound) g + random g = randomR (minBound,maxBound) g instance Random Bool where - randomR (a,b) g = + randomR (a,b) g = case (randomIvalInteger (bool2Int a, bool2Int b) g) of (x, g') -> (int2Bool x, g') where @@ -392,42 +392,42 @@ instance Random Bool where bool2Int False = 0 bool2Int True = 1 - int2Bool :: Int -> Bool - int2Bool 0 = False - int2Bool _ = True + int2Bool :: Int -> Bool + int2Bool 0 = False + int2Bool _ = True - random g = randomR (minBound,maxBound) g + random g = randomR (minBound,maxBound) g {-# INLINE randomRFloating #-} randomRFloating :: (Fractional a, Num a, Ord a, Random a, RandomGen g) => (a, a) -> g -> (a, g) -randomRFloating (l,h) g +randomRFloating (l,h) g | l>h = randomRFloating (h,l) g - | otherwise = let (coef,g') = random g in - (2.0 * (0.5*l + coef * (0.5*h - 0.5*l)), g') -- avoid overflow + | otherwise = let (coef,g') = random g in + (2.0 * (0.5*l + coef * (0.5*h - 0.5*l)), g') -- avoid overflow instance Random Double where randomR = randomRFloating - random rng = - case random rng of - (x,rng') -> + random rng = + case random rng of + (x,rng') -> -- We use 53 bits of randomness corresponding to the 53 bit significand: - ((fromIntegral (mask53 .&. (x::Int64)) :: Double) - / fromIntegral twoto53, rng') - where + ((fromIntegral (mask53 .&. (x::Int64)) :: Double) + / fromIntegral twoto53, rng') + where twoto53 = (2::Int64) ^ (53::Int64) mask53 = twoto53 - 1 - + instance Random Float where randomR = randomRFloating - random rng = - -- TODO: Faster to just use 'next' IF it generates enough bits of randomness. - case random rng of - (x,rng') -> + random rng = + -- TODO: Faster to just use 'next' IF it generates enough bits of randomness. + case random rng of + (x,rng') -> -- We use 24 bits of randomness corresponding to the 24 bit significand: - ((fromIntegral (mask24 .&. (x::Int32)) :: Float) - / fromIntegral twoto24, rng') - -- Note, encodeFloat is another option, but I'm not seeing slightly - -- worse performance with the following [2011.06.25]: + ((fromIntegral (mask24 .&. (x::Int32)) :: Float) + / fromIntegral twoto24, rng') + -- Note, encodeFloat is another option, but I'm not seeing slightly + -- worse performance with the following [2011.06.25]: -- (encodeFloat rand (-24), rng') where mask24 = twoto24 - 1 @@ -436,8 +436,8 @@ instance Random Float where -- CFloat/CDouble are basically the same as a Float/Double: instance Random CFloat where randomR = randomRFloating - random rng = case random rng of - (x,rng') -> (realToFrac (x::Float), rng') + random rng = case random rng of + (x,rng') -> (realToFrac (x::Float), rng') instance Random CDouble where randomR = randomRFloating @@ -445,8 +445,8 @@ instance Random CDouble where -- Presently, this is showing better performance than the Double instance: -- (And yet, if the Double instance uses randomFrac then its performance is much worse!) random = randomFrac - -- random rng = case random rng of - -- (x,rng') -> (realToFrac (x::Double), rng') + -- random rng = case random rng of + -- (x,rng') -> (realToFrac (x::Double), rng') mkStdRNG :: Integer -> IO StdGen mkStdRNG o = do @@ -463,7 +463,7 @@ randomIvalIntegral (l,h) = randomIvalInteger (toInteger l, toInteger h) {-# SPECIALIZE randomIvalInteger :: (Num a) => (Integer, Integer) -> StdGen -> (a, StdGen) #-} - + randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g) randomIvalInteger (l,h) rng | l > h = randomIvalInteger (h,l) rng @@ -482,7 +482,7 @@ randomIvalInteger (l,h) rng k = h - l + 1 magtgt = k * q - -- generate random values until we exceed the target magnitude + -- generate random values until we exceed the target magnitude f mag v g | mag >= magtgt = (v, g) | otherwise = v' `seq`f (mag*b) v' g' where (x,g') = next g @@ -494,18 +494,18 @@ randomFrac :: (RandomGen g, Fractional a) => g -> (a, g) randomFrac = randomIvalDouble (0::Double,1) realToFrac randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g) -randomIvalDouble (l,h) fromDouble rng +randomIvalDouble (l,h) fromDouble rng | l > h = randomIvalDouble (h,l) fromDouble rng - | otherwise = + | otherwise = case (randomIvalInteger (toInteger (minBound::Int32), toInteger (maxBound::Int32)) rng) of - (x, rng') -> - let - scaled_x = - fromDouble (0.5*l + 0.5*h) + -- previously (l+h)/2, overflowed + (x, rng') -> + let + scaled_x = + fromDouble (0.5*l + 0.5*h) + -- previously (l+h)/2, overflowed fromDouble ((0.5*h - 0.5*l) / (0.5 * realToFrac int32Count)) * -- avoid overflow - fromIntegral (x::Int32) - in - (scaled_x, rng') + fromIntegral (x::Int32) + in + (scaled_x, rng') int32Count :: Integer int32Count = toInteger (maxBound::Int32) - toInteger (minBound::Int32) + 1 -- GHC ticket #3982 @@ -516,16 +516,16 @@ 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' + 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) From 21520b265493a3eb64fd7769083c63a0a0f5d075 Mon Sep 17 00:00:00 2001 From: Herbert Valerio Riedel Date: Sun, 21 Dec 2014 10:45:03 +0100 Subject: [PATCH 02/31] Set GitHub's Git repo as primary repo --- random.cabal | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/random.cabal b/random.cabal index fd29840fb..a28063c90 100644 --- a/random.cabal +++ b/random.cabal @@ -40,7 +40,7 @@ Library source-repository head type: git - location: http://git.haskell.org/packages/random.git + location: https://github.com/haskell/random.git -- To run the Test-Suite: -- $ cabal configure --enable-tests From d9966ae732cb7e03dbd15b659fa852a7fad6692c Mon Sep 17 00:00:00 2001 From: Frank Cash Date: Thu, 12 Mar 2015 23:04:47 -0400 Subject: [PATCH 03/31] Update README.md Updates Haskell 98 link --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 9d5bb51b2..47b75896b 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ The API documentation can be found here: A module supplying this interface is required for Haskell 98 (but not Haskell 2010). An older [version] -(http://www.haskell.org/ghc/docs/latest/html/libraries/haskell98/Random.html) +(https://downloads.haskell.org/~ghc/latest/docs/html/libraries/haskell98-2.0.0.3/Random.html) of this library is included with GHC in the haskell98 package. This newer version, with compatible api, is included in the [Haskell Platform] (http://www.haskell.org/platform/contents.html). From 5af0c12ed4e844400cf89c8f139b88a3e6f02067 Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Sat, 4 Apr 2015 17:29:03 +0100 Subject: [PATCH 04/31] Add some tests demonstrating https://github.com/haskell/random/issues/25. --- Benchmark/RandomnessTest2.hs | 42 +++++++++++++++++++++++++++++++++++ Benchmark/RandomnessTest3.hs | 43 ++++++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+) create mode 100644 Benchmark/RandomnessTest2.hs create mode 100644 Benchmark/RandomnessTest3.hs diff --git a/Benchmark/RandomnessTest2.hs b/Benchmark/RandomnessTest2.hs new file mode 100644 index 000000000..14bc0d773 --- /dev/null +++ b/Benchmark/RandomnessTest2.hs @@ -0,0 +1,42 @@ +module Main where + +import Control.Monad +import System.Random + +-- splitN :: (RandomGen g) => Int -> g -> ([g], g) +splitN 0 g = ([], g) +splitN n g = (g1:l, g') + where + (l, g') = splitN (n-1) g2 + (g1, g2) = split g + +-- The funny splitting operation. +split' :: (RandomGen g) => g -> (g, g) +split' g = (g12, g21) + where + (g1, g2) = split g + (_, g12) = split g1 + (g21, _) = split g2 + +-- This test checks if generators created by calling split 2 times are independent. +-- It generates pairs of integers from 0 to n-1, using split' to +-- generate both numbers using one seed. Then it counts how often the +-- two numbers are equal. +test :: (RandomGen g) => Int -> Int -> g -> Int +test numTests n g = equals + where + (gs, _) = splitN numTests g + equals = count id $ map single gs + count p l = length $ filter p l + single g' = (fst $ randomR (0, n-1) g1) == (fst $ randomR (0, n-1) g2) + where + (g1, g2) = split' g' + +main = do + let g = mkStdGen 42 + forM_ [2..15] $ \i -> do + let actual = test (i * 1000) i g + putStrLn $ "Generated " ++ show (i * 1000) + ++ " pairs of numbers from 0 to " ++ show (i - 1) + ++ " -- " ++ show actual ++ " pairs contained equal numbers " + ++ "and we expected about 1000." diff --git a/Benchmark/RandomnessTest3.hs b/Benchmark/RandomnessTest3.hs new file mode 100644 index 000000000..b943f59d3 --- /dev/null +++ b/Benchmark/RandomnessTest3.hs @@ -0,0 +1,43 @@ +module Main where + +import Control.Monad +import System.Random.TF.Gen +import System.Random.TF.Instances + +splitN :: (RandomGen g) => Int -> g -> ([g], g) +splitN 0 g = ([], g) +splitN n g = (g1:l, g') + where + (l, g') = splitN (n-1) g2 + (g1, g2) = split g + +-- The funny splitting operation. +split' :: (RandomGen g) => g -> (g, g) +split' g = (g12, g21) + where + (g1, g2) = split g + (_, g12) = split g1 + (g21, _) = split g2 + +-- This test checks if generators created by calling split 2 times are independent. +-- It generates pairs of integers from 0 to n-1, using split' to +-- generate both numbers using one seed. Then it counts how often the +-- two numbers are equal. +test :: (RandomGen g) => Int -> Int -> g -> Int +test numTests n g = equals + where + (gs, _) = splitN numTests g + equals = count id $ map single gs + count p l = length $ filter p l + single g' = (fst $ randomR (0, n-1) g1) == (fst $ randomR (0, n-1) g2) + where + (g1, g2) = split' g' + +main = do + let g = seedTFGen (42, 42, 42, 42) + forM_ [2..15] $ \i -> do + let actual = test (i * 1000) i g + putStrLn $ "Generated " ++ show (i * 1000) + ++ " pairs of numbers from 0 to " ++ show (i - 1) + ++ " -- " ++ show actual ++ " pairs contained equal numbers " + ++ "and we expected about 1000." From d88796ea55a53bafe6b1bbbd95a3555f4e122952 Mon Sep 17 00:00:00 2001 From: Dominic Steinitz Date: Sat, 4 Apr 2015 18:25:33 +0100 Subject: [PATCH 05/31] Make benchmark build and add benchmark for tf-random. --- Benchmark/Makefile | 3 +-- Benchmark/SimpleRNGBench.hs | 13 ++++++++++--- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/Benchmark/Makefile b/Benchmark/Makefile index 8a84e6479..69d2d37b5 100644 --- a/Benchmark/Makefile +++ b/Benchmark/Makefile @@ -1,7 +1,6 @@ -#OPTS= -O2 -Wall -XCPP -OPTS= -O2 -Wall -XCPP -Werror +OPTS= -O2 -Wall -XCPP all: lib bench diff --git a/Benchmark/SimpleRNGBench.hs b/Benchmark/SimpleRNGBench.hs index c25b75d47..04e2f7715 100644 --- a/Benchmark/SimpleRNGBench.hs +++ b/Benchmark/SimpleRNGBench.hs @@ -33,6 +33,7 @@ import Prelude hiding (last,sum) import BinSearch #ifdef TEST_COMPETITORS +import System.Random.TF import System.Random.Mersenne.Pure64 import System.Random.MWC import Control.Monad.Primitive @@ -44,11 +45,11 @@ import GHC.IO -- Miscellaneous helpers: -- Readable large integer printing: -commaint :: Integral a => a -> String +commaint :: (Show a, Integral a) => a -> String commaint n = reverse $ concat $ intersperse "," $ - chunk 3 $ + chunksOf 3 $ reverse (show n) padleft :: Int -> String -> String @@ -281,7 +282,13 @@ main = do timeit th freq "System.Random.Mersenne.Pure64 Ints" gen_mt randInt2 timeit th freq "System.Random.Mersenne.Pure64 Floats" gen_mt randFloat2 --- gen_mwc <- create + let gen_tf = seedTFGen (0,1,2,3) + randInt4 = random :: RandomGen g => g -> (Int,g) + randFloat4 = random :: RandomGen g => g -> (Float,g) + timeit th freq "System.Random.TF next" gen_tf next + timeit th freq "System.Random.TF Ints" gen_tf randInt4 + timeit th freq "System.Random.TF Floats" gen_tf randFloat4 + withSystemRandom $ \ gen_mwc -> do let randInt3 = random :: RandomGen g => g -> (Int,g) randFloat3 = random :: RandomGen g => g -> (Float,g) From 357090b391d2bc2ada54f4f2c8592f7f854f60cb Mon Sep 17 00:00:00 2001 From: Thomas Miedema Date: Wed, 20 Jan 2016 16:35:06 +0100 Subject: [PATCH 06/31] GHC testsuite: require random to be installed --- tests/all.T | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/all.T b/tests/all.T index f1675ed5c..f85ec767a 100644 --- a/tests/all.T +++ b/tests/all.T @@ -1,3 +1,4 @@ +setTestOpts(reqlib('random')) test('rangeTest', normal, From bc34220740c467f0154d3da0eefc2b3a2db6ae96 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sun, 15 Jan 2017 14:43:11 -0500 Subject: [PATCH 07/31] version 2.0 of random begins --- .darcs-boring | 5 - .gitignore | 12 - .travis.yml | 5 - Benchmark/BinSearch.hs | 142 -------- Benchmark/Makefile | 23 -- Benchmark/RandomnessTest2.hs | 42 --- Benchmark/RandomnessTest3.hs | 43 --- Benchmark/SimpleRNGBench.hs | 329 ------------------- CHANGELOG.md | 27 +- DEVLOG.md | 196 ----------- LICENSE | 83 ++--- README.md | 18 -- Setup.hs | 4 - System/Random.hs | 609 ----------------------------------- prologue.txt | 1 - random.cabal | 140 ++++---- tests/Makefile | 14 - tests/T7936.hs | 14 - tests/TestRandomIOs.hs | 20 -- tests/TestRandomRs.hs | 22 -- tests/all.T | 11 - tests/random1283.hs | 40 --- tests/rangeTest.hs | 135 -------- tests/rangeTest.stdout | 67 ---- tests/slowness.hs | 55 ---- 25 files changed, 96 insertions(+), 1961 deletions(-) delete mode 100644 .darcs-boring delete mode 100644 .gitignore delete mode 100644 .travis.yml delete mode 100644 Benchmark/BinSearch.hs delete mode 100644 Benchmark/Makefile delete mode 100644 Benchmark/RandomnessTest2.hs delete mode 100644 Benchmark/RandomnessTest3.hs delete mode 100644 Benchmark/SimpleRNGBench.hs delete mode 100644 DEVLOG.md delete mode 100644 README.md delete mode 100644 System/Random.hs delete mode 100644 prologue.txt delete mode 100644 tests/Makefile delete mode 100644 tests/T7936.hs delete mode 100644 tests/TestRandomIOs.hs delete mode 100644 tests/TestRandomRs.hs delete mode 100644 tests/all.T delete mode 100644 tests/random1283.hs delete mode 100644 tests/rangeTest.hs delete mode 100644 tests/rangeTest.stdout delete mode 100644 tests/slowness.hs diff --git a/.darcs-boring b/.darcs-boring deleted file mode 100644 index 6682606a0..000000000 --- a/.darcs-boring +++ /dev/null @@ -1,5 +0,0 @@ -^dist(/|$) -^setup(/|$) -^GNUmakefile$ -^Makefile.local$ -^.depend(.bak)?$ diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 41c6d8c61..000000000 --- a/.gitignore +++ /dev/null @@ -1,12 +0,0 @@ -*~ - -Thumbs.db -.DS_Store - -GNUmakefile -dist-install/ -ghc.mk - -dist -.cabal-sandbox -cabal.sandbox.config diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 6a03bcb12..000000000 --- a/.travis.yml +++ /dev/null @@ -1,5 +0,0 @@ -language: haskell -ghc: - - 7.4 - - 7.6 - - 7.8 diff --git a/Benchmark/BinSearch.hs b/Benchmark/BinSearch.hs deleted file mode 100644 index f61164855..000000000 --- a/Benchmark/BinSearch.hs +++ /dev/null @@ -1,142 +0,0 @@ - -{- - Binary search over benchmark input sizes. - - There are many good ways to measure the time it takes to perform a - certain computation on a certain input. However, frequently, it's - challenging to pick the right input size for all platforms and all - compilataion modes. - - Sometimes for linear-complexity benchmarks it is better to measure - /throughput/, i.e. elements processed per second. That is, fixing - the time of execution and measuring the amount of work done (rather - than the reverse). This library provides a simple way to search for - an appropriate input size that results in the desired execution time. - - An alternative approach is to kill the computation after a certain - amount of time and observe how much work it has completed. - -} -module BinSearch - ( - binSearch - ) -where - -import Control.Monad -import Data.Time.Clock -- Not in 6.10 -import Data.List -import System.IO -import Prelude hiding (min,max,log) - - - --- | Binary search for the number of inputs to a computation that --- results in a specified amount of execution time in seconds. For example: --- --- > binSearch verbose N (min,max) kernel --- --- ... will find the right input size that results in a time --- between min and max, then it will then run for N trials and --- return the median (input,time-in-seconds) pair. -binSearch :: Bool -> Integer -> (Double,Double) -> (Integer -> IO ()) -> IO (Integer, Double) -binSearch verbose trials (min,max) kernel = - do - when(verbose)$ putStrLn$ "[binsearch] Binary search for input size resulting in time in range "++ show (min,max) - - let desired_exec_length = 1.0 - good_trial t = (toRational t <= toRational max) && (toRational t >= toRational min) - - -- At some point we must give up... - loop n | n > ((2::Integer) ^ (100::Integer)) = error "ERROR binSearch: This function doesn't seem to scale in proportion to its last argument." - - -- Not allowed to have "0" size input, bump it back to one: - loop 0 = loop 1 - - loop n = - do - when(verbose)$ putStr$ "[binsearch:"++ show n ++ "] " - time <- timeit$ kernel n - when(verbose)$ putStrLn$ "Time consumed: "++ show time - let rate = fromIntegral n / time - - -- [2010.06.09] Introducing a small fudge factor to help our guess get over the line: - let initial_fudge_factor = 1.10 - fudge_factor = 1.01 -- Even in the steady state we fudge a little - guess = desired_exec_length * rate - - -- TODO: We should keep more history here so that we don't re-explore input space we have already explored. - -- This is a balancing act because of randomness in execution time. - - if good_trial time - then do - when(verbose)$ putStrLn$ "[binsearch] Time in range. LOCKING input size and performing remaining trials." - print_trial 1 n time - lockin (trials-1) n [time] - - -- Here we're still in the doubling phase: - else if time < 0.100 - then loop (2*n) - - else do when(verbose)$ - putStrLn$ "[binsearch] Estimated rate to be " - ++show (round rate::Integer)++" per second. Trying to scale up..." - - -- Here we've exited the doubling phase, but we're making our first guess as to how big a real execution should be: - if time > 0.100 && time < 0.33 * desired_exec_length - then do when(verbose)$ putStrLn$ "[binsearch] (Fudging first guess a little bit extra)" - loop (round$ guess * initial_fudge_factor) - else loop (round$ guess * fudge_factor) - - -- Termination condition: Done with all trials. - lockin 0 n log = do when(verbose)$ putStrLn$ "[binsearch] Time-per-unit for all trials: "++ - (concat $ intersperse " " (map (show . (/ toDouble n) . toDouble) $ sort log)) - return (n, log !! ((length log) `quot` 2)) -- Take the median - - lockin trials_left n log = - do when(verbose)$ putStrLn$ "[binsearch]------------------------------------------------------------" - time <- timeit$ kernel n - -- hFlush stdout - print_trial (trials - trials_left +1 ) n time - -- when(verbose)$ hFlush stdout - lockin (trials_left - 1) n (time : log) - - print_trial :: Integer -> Integer -> NominalDiffTime -> IO () - print_trial trialnum n time = - let rate = fromIntegral n / time - timeperunit = time / fromIntegral n - in - when(verbose)$ putStrLn$ "[binsearch] TRIAL: "++show trialnum ++ - " secPerUnit: "++ showTime timeperunit ++ - " ratePerSec: "++ show (rate) ++ - " seconds: "++showTime time - - - - (n,t) <- loop 1 - return (n, fromRational$ toRational t) - -showTime :: NominalDiffTime -> String -showTime t = show ((fromRational $ toRational t) :: Double) - -toDouble :: Real a => a -> Double -toDouble = fromRational . toRational - - --- Could use cycle counters here.... but the point of this is to time --- things on the order of a second. -timeit :: IO () -> IO NominalDiffTime -timeit io = - do strt <- getCurrentTime - io - end <- getCurrentTime - return (diffUTCTime end strt) -{- -test :: IO (Integer,Double) -test = - binSearch True 3 (1.0, 1.05) - (\n -> - do v <- newIORef 0 - forM_ [1..n] $ \i -> do - old <- readIORef v - writeIORef v (old+i)) --} diff --git a/Benchmark/Makefile b/Benchmark/Makefile deleted file mode 100644 index 69d2d37b5..000000000 --- a/Benchmark/Makefile +++ /dev/null @@ -1,23 +0,0 @@ - - -OPTS= -O2 -Wall -XCPP - -all: lib bench - -lib: - (cd .. && ghc $(OPTS) --make System/Random.hs) - -bench: - ghc $(OPTS) -i.. --make SimpleRNGBench.hs - -# When benchmarking against other packages installed via cabal it is -# necessary to IGNORE the System/Random.hs files in the current directory. -# (Otherwise instances of RandomGen are confused.) -bench2: - ghc $(OPTS) -DTEST_COMPETITORS --make SimpleRNGBench.hs - -clean: - rm -f *.o *.hi SimpleRNGBench -# cabal clean -# (cd Benchmark/ && rm -f *.o *.hi SimpleRNGBench) -# (cd System/ && rm -f *.o *.hi SimpleRNGBench) diff --git a/Benchmark/RandomnessTest2.hs b/Benchmark/RandomnessTest2.hs deleted file mode 100644 index 14bc0d773..000000000 --- a/Benchmark/RandomnessTest2.hs +++ /dev/null @@ -1,42 +0,0 @@ -module Main where - -import Control.Monad -import System.Random - --- splitN :: (RandomGen g) => Int -> g -> ([g], g) -splitN 0 g = ([], g) -splitN n g = (g1:l, g') - where - (l, g') = splitN (n-1) g2 - (g1, g2) = split g - --- The funny splitting operation. -split' :: (RandomGen g) => g -> (g, g) -split' g = (g12, g21) - where - (g1, g2) = split g - (_, g12) = split g1 - (g21, _) = split g2 - --- This test checks if generators created by calling split 2 times are independent. --- It generates pairs of integers from 0 to n-1, using split' to --- generate both numbers using one seed. Then it counts how often the --- two numbers are equal. -test :: (RandomGen g) => Int -> Int -> g -> Int -test numTests n g = equals - where - (gs, _) = splitN numTests g - equals = count id $ map single gs - count p l = length $ filter p l - single g' = (fst $ randomR (0, n-1) g1) == (fst $ randomR (0, n-1) g2) - where - (g1, g2) = split' g' - -main = do - let g = mkStdGen 42 - forM_ [2..15] $ \i -> do - let actual = test (i * 1000) i g - putStrLn $ "Generated " ++ show (i * 1000) - ++ " pairs of numbers from 0 to " ++ show (i - 1) - ++ " -- " ++ show actual ++ " pairs contained equal numbers " - ++ "and we expected about 1000." diff --git a/Benchmark/RandomnessTest3.hs b/Benchmark/RandomnessTest3.hs deleted file mode 100644 index b943f59d3..000000000 --- a/Benchmark/RandomnessTest3.hs +++ /dev/null @@ -1,43 +0,0 @@ -module Main where - -import Control.Monad -import System.Random.TF.Gen -import System.Random.TF.Instances - -splitN :: (RandomGen g) => Int -> g -> ([g], g) -splitN 0 g = ([], g) -splitN n g = (g1:l, g') - where - (l, g') = splitN (n-1) g2 - (g1, g2) = split g - --- The funny splitting operation. -split' :: (RandomGen g) => g -> (g, g) -split' g = (g12, g21) - where - (g1, g2) = split g - (_, g12) = split g1 - (g21, _) = split g2 - --- This test checks if generators created by calling split 2 times are independent. --- It generates pairs of integers from 0 to n-1, using split' to --- generate both numbers using one seed. Then it counts how often the --- two numbers are equal. -test :: (RandomGen g) => Int -> Int -> g -> Int -test numTests n g = equals - where - (gs, _) = splitN numTests g - equals = count id $ map single gs - count p l = length $ filter p l - single g' = (fst $ randomR (0, n-1) g1) == (fst $ randomR (0, n-1) g2) - where - (g1, g2) = split' g' - -main = do - let g = seedTFGen (42, 42, 42, 42) - forM_ [2..15] $ \i -> do - let actual = test (i * 1000) i g - putStrLn $ "Generated " ++ show (i * 1000) - ++ " pairs of numbers from 0 to " ++ show (i - 1) - ++ " -- " ++ show actual ++ " pairs contained equal numbers " - ++ "and we expected about 1000." diff --git a/Benchmark/SimpleRNGBench.hs b/Benchmark/SimpleRNGBench.hs deleted file mode 100644 index 04e2f7715..000000000 --- a/Benchmark/SimpleRNGBench.hs +++ /dev/null @@ -1,329 +0,0 @@ -{-# LANGUAGE BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-} -{-# OPTIONS_GHC -fwarn-unused-imports #-} - --- | A simple script to do some very basic timing of the RNGs. - -module Main where - -import System.Exit (exitSuccess, exitFailure) -import System.Environment -import System.Random -import System.CPUTime (getCPUTime) -import System.CPUTime.Rdtsc -import System.Console.GetOpt - -import GHC.Conc -import Control.Concurrent -import Control.Monad -import Control.Exception - -import Data.IORef -import Data.Word -import Data.List hiding (last,sum) -import Data.Int -import Data.List.Split hiding (split) -import Text.Printf - -import Foreign.Ptr -import Foreign.C.Types -import Foreign.ForeignPtr -import Foreign.Storable (peek,poke) - -import Prelude hiding (last,sum) -import BinSearch - -#ifdef TEST_COMPETITORS -import System.Random.TF -import System.Random.Mersenne.Pure64 -import System.Random.MWC -import Control.Monad.Primitive --- import System.IO.Unsafe -import GHC.IO -#endif - ----------------------------------------------------------------------------------------------------- --- Miscellaneous helpers: - --- Readable large integer printing: -commaint :: (Show a, Integral a) => a -> String -commaint n = - reverse $ concat $ - intersperse "," $ - chunksOf 3 $ - reverse (show n) - -padleft :: Int -> String -> String -padleft n str | length str >= n = str -padleft n str | otherwise = take (n - length str) (repeat ' ') ++ str - -padright :: Int -> String -> String -padright n str | length str >= n = str -padright n str | otherwise = str ++ take (n - length str) (repeat ' ') - -fmt_num :: (RealFrac a, PrintfArg a) => a -> String -fmt_num n = if n < 100 - then printf "%.2f" n - else commaint (round n :: Integer) - - --- Measure clock frequency, spinning rather than sleeping to try to --- stay on the same core. -measureFreq :: IO Int64 -measureFreq = do - let second = 1000 * 1000 * 1000 * 1000 -- picoseconds are annoying - t1 <- rdtsc - start <- getCPUTime - let loop !n !last = - do t2 <- rdtsc - when (t2 < last) $ - putStrLn$ "COUNTERS WRAPPED "++ show (last,t2) - cput <- getCPUTime - if (cput - start < second) - then loop (n+1) t2 - else return (n,t2) - (n,t2) <- loop 0 t1 - putStrLn$ " Approx getCPUTime calls per second: "++ commaint (n::Int64) - when (t2 < t1) $ - putStrLn$ "WARNING: rdtsc not monotonically increasing, first "++show t1++" then "++show t2++" on the same OS thread" - - return$ fromIntegral (t2 - t1) - ----------------------------------------------------------------------------------------------------- - --- Test overheads without actually generating any random numbers: -data NoopRNG = NoopRNG -instance RandomGen NoopRNG where - next g = (0,g) -#ifdef ENABLE_SPLITTABLEGEN - genRange _ = (0,0) -instance SplittableGen NoopRNG where -#endif - split g = (g,g) - --- An RNG generating only 0 or 1: -data BinRNG = BinRNG StdGen -instance RandomGen BinRNG where - next (BinRNG g) = (x `mod` 2, BinRNG g') - where (x,g') = next g -#ifdef ENABLE_SPLITTABLEGEN - genRange _ = (0,1) -instance SplittableGen BinRNG where -#endif - split (BinRNG g) = (BinRNG g1, BinRNG g2) - where (g1,g2) = split g - - - -#ifdef TEST_COMPETITORS -data MWCRNG = MWCRNG (Gen (PrimState IO)) --- data MWCRNG = MWCRNG GenIO -instance RandomGen MWCRNG where - -- For testing purposes we hack this to be non-monadic: --- next g@(MWCRNG gen) = unsafePerformIO $ - next g@(MWCRNG gen) = unsafeDupablePerformIO $ - do v <- uniform gen - return (v, g) -#endif - ----------------------------------------------------------------------------------------------------- --- Drivers to get random numbers repeatedly. - -type Kern = Int -> Ptr Int -> IO () - --- [2011.01.28] Changing this to take "count" and "accumulator ptr" as arguments: --- foreign import ccall "cbits/c_test.c" blast_rands :: Kern --- foreign import ccall "cbits/c_test.c" store_loop :: Kern --- foreign import ccall unsafe "stdlib.hs" rand :: IO Int - -{-# INLINE timeit #-} -timeit :: (Random a, RandomGen g) => - Int -> Int64 -> String -> g -> (g -> (a,g)) -> IO () -timeit numthreads freq msg gen nxt = - do - counters <- forM [1..numthreads] (const$ newIORef (1::Int64)) - tids <- forM counters $ \counter -> - forkIO $ infloop counter (nxt gen) - threadDelay (1000*1000) -- One second - mapM_ killThread tids - - finals <- mapM readIORef counters - let mean :: Double = fromIntegral (foldl1 (+) finals) / fromIntegral numthreads - cycles_per :: Double = fromIntegral freq / mean - printResult (round mean :: Int64) msg cycles_per - - where - infloop !counter !(!_,!g) = - do incr counter - infloop counter (nxt g) - - incr !counter = - do -- modifyIORef counter (+1) -- Not strict enough! - c <- readIORef counter - let c' = c+1 - _ <- evaluate c' - writeIORef counter c' - - --- This function times an IO function on one or more threads. Rather --- than running a fixed number of iterations, it uses a binary search --- to find out how many iterations can be completed in a second. -timeit_foreign :: Int -> Int64 -> String -> (Int -> Ptr Int -> IO ()) -> IO Int64 -timeit_foreign numthreads freq msg ffn = do - ptr :: ForeignPtr Int <- mallocForeignPtr - - let kern = if numthreads == 1 - then ffn - else replicate_kernel numthreads ffn - wrapped n = withForeignPtr ptr (kern$ fromIntegral n) - (n,t) <- binSearch False 1 (1.0, 1.05) wrapped - - let total_per_second = round $ fromIntegral n * (1 / t) - cycles_per = fromIntegral freq * t / fromIntegral n - printResult total_per_second msg cycles_per - return total_per_second - - where - -- This lifts a C kernel to operate simultaneously on N threads. - replicate_kernel :: Int -> Kern -> Kern - replicate_kernel nthreads kern n ptr = do - ptrs <- forM [1..nthreads] - (const mallocForeignPtr) - tmpchan <- newChan - -- let childwork = ceiling$ fromIntegral n / fromIntegral nthreads - let childwork = n -- Keep it the same.. interested in per-thread throughput. - -- Fork/join pattern: - _ <- forM ptrs $ \pt -> forkIO $ - withForeignPtr pt $ \p -> do - kern (fromIntegral childwork) p - result <- peek p - writeChan tmpchan result - - results <- forM [1..nthreads] $ \_ -> - readChan tmpchan - -- Meaningless semantics here... sum the child ptrs and write to the input one: - poke ptr (foldl1 (+) results) - return () - - -printResult :: Int64 -> String -> Double -> IO () -printResult total msg cycles_per = - putStrLn$ " "++ padleft 11 (commaint total) ++" randoms generated "++ padright 27 ("["++msg++"]") ++" ~ " - ++ fmt_num cycles_per ++" cycles/int" - ----------------------------------------------------------------------------------------------------- --- Main Script - -data Flag = NoC | Help - deriving (Show, Eq) - -options :: [OptDescr Flag] -options = - [ Option ['h'] ["help"] (NoArg Help) "print program help" - , Option [] ["noC"] (NoArg NoC) "omit C benchmarks, haskell only" - ] - -main :: IO () -main = do - argv <- getArgs - let (opts,_,other) = getOpt Permute options argv - - when (not$ null other) $ do - putStrLn$ "ERROR: Unrecognized options: " - mapM_ putStr other - exitFailure - - when (Help `elem` opts) $ do - putStr$ usageInfo "Benchmark random number generation" options - exitSuccess - - putStrLn$ "\nHow many random numbers can we generate in a second on one thread?" - - t1 <- rdtsc - t2 <- rdtsc - putStrLn (" Cost of rdtsc (ffi call): " ++ show (t2 - t1)) - - freq <- measureFreq - putStrLn$ " Approx clock frequency: " ++ commaint freq - - let - randInt = random :: RandomGen g => g -> (Int,g) - randWord16 = random :: RandomGen g => g -> (Word16,g) - randFloat = random :: RandomGen g => g -> (Float,g) - randCFloat = random :: RandomGen g => g -> (CFloat,g) - randDouble = random :: RandomGen g => g -> (Double,g) - randCDouble = random :: RandomGen g => g -> (CDouble,g) - randInteger = random :: RandomGen g => g -> (Integer,g) - randBool = random :: RandomGen g => g -> (Bool,g) - randChar = random :: RandomGen g => g -> (Char,g) - - gen = mkStdGen 23852358661234 - gamut th = do - putStrLn$ " First, timing System.Random.next:" - timeit th freq "constant zero gen" NoopRNG next - timeit th freq "System.Random stdGen/next" gen next - - putStrLn$ "\n Second, timing System.Random.random at different types:" - timeit th freq "System.Random Ints" gen randInt - timeit th freq "System.Random Word16" gen randWord16 - timeit th freq "System.Random Floats" gen randFloat - timeit th freq "System.Random CFloats" gen randCFloat - timeit th freq "System.Random Doubles" gen randDouble - timeit th freq "System.Random CDoubles" gen randCDouble - timeit th freq "System.Random Integers" gen randInteger - timeit th freq "System.Random Bools" gen randBool - timeit th freq "System.Random Chars" gen randChar - -#ifdef TEST_COMPETITORS - putStrLn$ "\n Next test other RNG packages on Hackage:" - let gen_mt = pureMT 39852 - randInt2 = random :: RandomGen g => g -> (Int,g) - randFloat2 = random :: RandomGen g => g -> (Float,g) - timeit th freq "System.Random.Mersenne.Pure64 next" gen_mt next - timeit th freq "System.Random.Mersenne.Pure64 Ints" gen_mt randInt2 - timeit th freq "System.Random.Mersenne.Pure64 Floats" gen_mt randFloat2 - - let gen_tf = seedTFGen (0,1,2,3) - randInt4 = random :: RandomGen g => g -> (Int,g) - randFloat4 = random :: RandomGen g => g -> (Float,g) - timeit th freq "System.Random.TF next" gen_tf next - timeit th freq "System.Random.TF Ints" gen_tf randInt4 - timeit th freq "System.Random.TF Floats" gen_tf randFloat4 - - withSystemRandom $ \ gen_mwc -> do - let randInt3 = random :: RandomGen g => g -> (Int,g) - randFloat3 = random :: RandomGen g => g -> (Float,g) - - timeit th freq "System.Random.MWC next" (MWCRNG gen_mwc) next - timeit th freq "System.Random.MWC Ints" (MWCRNG gen_mwc) randInt3 - timeit th freq "System.Random.MWC Floats" (MWCRNG gen_mwc) randFloat3 - -#endif - - putStrLn$ "\n Next timing range-restricted System.Random.randomR:" - timeit th freq "System.Random Ints" gen (randomR (-100, 100::Int)) - timeit th freq "System.Random Word16s" gen (randomR (-100, 100::Word16)) - timeit th freq "System.Random Floats" gen (randomR (-100, 100::Float)) - timeit th freq "System.Random CFloats" gen (randomR (-100, 100::CFloat)) - timeit th freq "System.Random Doubles" gen (randomR (-100, 100::Double)) - timeit th freq "System.Random CDoubles" gen (randomR (-100, 100::CDouble)) - timeit th freq "System.Random Integers" gen (randomR (-100, 100::Integer)) - timeit th freq "System.Random Bools" gen (randomR (False, True::Bool)) - timeit th freq "System.Random Chars" gen (randomR ('a', 'z')) - timeit th freq "System.Random BIG Integers" gen (randomR (0, (2::Integer) ^ (5000::Int))) - - -- when (not$ NoC `elem` opts) $ do - -- putStrLn$ " Comparison to C's rand():" - -- timeit_foreign th freq "ptr store in C loop" store_loop - -- timeit_foreign th freq "rand/store in C loop" blast_rands - -- timeit_foreign th freq "rand in Haskell loop" (\n ptr -> forM_ [1..n]$ \_ -> rand ) - -- timeit_foreign th freq "rand/store in Haskell loop" (\n ptr -> forM_ [1..n]$ \_ -> do n <- rand; poke ptr n ) - -- return () - - -- Test with 1 thread and numCapabilities threads: - gamut 1 - when (numCapabilities > 1) $ do - putStrLn$ "\nNow "++ show numCapabilities ++" threads, reporting mean randoms-per-second-per-thread:" - gamut numCapabilities - return () - - putStrLn$ "Finished." diff --git a/CHANGELOG.md b/CHANGELOG.md index 15c882af9..71040582c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,26 +1,5 @@ -# 1.1 - * breaking change to `randomIValInteger` to improve RNG quality and performance - see https://github.com/haskell/random/pull/4 and - ghc https://ghc.haskell.org/trac/ghc/ticket/8898 - * correct documentation about generated range of Int32 sized values of type Int - https://github.com/haskell/random/pull/7 - * fix memory leaks by using strict fields and strict atomicModifyIORef' - https://github.com/haskell/random/pull/8 - related to ghc trac tickets #7936 and #4218 - * support for base < 4.6 (which doesnt provide strict atomicModifyIORef') - and integrating Travis CI support. - https://github.com/haskell/random/pull/12 - * fix C type in test suite https://github.com/haskell/random/pull/9 +# Revision history for random -# 1.0.1.1 -bump for overflow bug fixes - -# 1.0.1.2 -bump for ticket 8704, build fusion - -# 1.0.1.0 -bump for bug fixes, - -# 1.0.0.4 -bumped version for float/double range bugfix +## 2.0.0.0 -- YYYY-mm-dd +* First version. Released on an unsuspecting world. diff --git a/DEVLOG.md b/DEVLOG.md deleted file mode 100644 index 6e0b28dc4..000000000 --- a/DEVLOG.md +++ /dev/null @@ -1,196 +0,0 @@ - - -DEVLOG: A collection of notes accumulated during development. -============================================================= - - -[2011.06.24] (transient) Regression in stdGen performance. ----------------------------------------------------------- - -I just added a simple benchmark to make sure that whatever fix I -introduce for trac ticket #5133 does not regress performance. Yet in -doing so I discovered that I'm getting much worse performance out of -rev 130e421e912d than I'm seeing in my installed random-1.0.0.3 package. - -Current version: - How many random numbers can we generate in a second on one thread? - Cost of rdtsc (ffi call): 100 - Approx getCPUTime calls per second: 234,553 - Approx clock frequency: 3,335,220,196 - First, timing with System.Random interface: - 68,550,189 random ints generated [constant zero gen] ~ 48.65 cycles/int - 900,889 random ints generated [System.Random stdGen] ~ 3,702 cycles/int - -random-1.0.0.3 version: - How many random numbers can we generate in a second on one thread? - Cost of rdtsc (ffi call): 75 - Approx getCPUTime calls per second: 215,332 - Approx clock frequency: 3,334,964,738 - First, timing with System.Random interface: - 71,683,748 random ints generated [constant zero gen] ~ 46.52 cycles/int - 13,609,559 random ints generated [System.Random stdGen] ~ 245 cycles/int - -A >13X difference!! -Both are compiled with the same options. The only difference is which -System.Random is used. - -When did the regression occur? - - * e059ed15172585310f9c -- 10/13/2010 perf still good - * 6c43f80f48178ac617 -- SplittableGen introduced, still good perf - * 130e421e912d394653a4 -- most recent, bad performance - -Ok... this is very odd. It was a heisenbug becuase it's disappeared -now! I'll leave this note here to help remember to look for it in the -future. - -Ryan - - -[2011.06.24] Timing non-int types ---------------------------------- - -The results are highly uneven: - - Cost of rdtsc (ffi call): 84 - Approx getCPUTime calls per second: 220,674 - Approx clock frequency: 3,336,127,273 - First, timing with System.Random interface: - 112,976,933 randoms generated [constant zero gen] ~ 29.53 cycles/int - 14,415,176 randoms generated [System.Random stdGen] ~ 231 cycles/int - 70,751 randoms generated [System.Random Floats] ~ 47,153 cycles/int - 70,685 randoms generated [System.Random CFloats] ~ 47,197 cycles/int - 2,511,635 randoms generated [System.Random Doubles] ~ 1,328 cycles/int - 70,494 randoms generated [System.Random CDoubles] ~ 47,325 cycles/int - 858,012 randoms generated [System.Random Integers] ~ 3,888 cycles/int - 4,756,213 randoms generated [System.Random Bools] ~ 701 cycles/int - -As you can see, all the types that use the generic randomIvalFrac / -randomFrac definitions perform badly. What's more, the above results -INCLUDE an attempt to inline: - - {-# INLINE randomIvalFrac #-} - {-# INLINE randomFrac #-} - {-# INLINE randomIvalDouble #-} - -After reimplementing random/Float these are the new results: - - Cost of rdtsc (ffi call): 100 - Approx getCPUTime calls per second: 200,582 - Approx clock frequency: 3,334,891,942 - First, timing with System.Random interface: - 105,266,949 randoms generated [constant zero gen] ~ 31.68 cycles/int - 13,593,392 randoms generated [System.Random stdGen] ~ 245 cycles/int - 10,962,597 randoms generated [System.Random Floats] ~ 304 cycles/int - 11,926,573 randoms generated [System.Random CFloats] ~ 280 cycles/int - 2,421,520 randoms generated [System.Random Doubles] ~ 1,377 cycles/int - 2,535,087 randoms generated [System.Random CDoubles] ~ 1,315 cycles/int - 856,276 randoms generated [System.Random Integers] ~ 3,895 cycles/int - 4,976,373 randoms generated [System.Random Bools] ~ 670 cycles/int - -(But I still need to propagate these changes throughout all types / API calls.) - - - -[2011.06.28] Integer Generation via random and randomR -------------------------------------------------------- - -Back on the master branch I notice that while randomIvalInteger does -well for small ranges, it's advantage doesn't scale to larger ranges: - - range (-100,100): - 5,105,290 randoms generated [System.Random Integers] ~ 653 cycles/int - - range (0,2^5000): - 8,969 randoms generated [System.Random BIG Integers] ~ 371,848 cycles/int - - - -[2011.08.25] Validating release version 1.0.1.0 rev 40bbfd2867 --------------------------------------------------------------- - -This is a bugfix release without SplittableGen. It passed (cd tests; -make test) on my Mac Os 10.6 machine. - -I ran GHC validate using the following fingerprint - - .|c5056b932a06b4adce5167a5cb69f1f0768d28ec - ghc-tarballs|e7b7b152083f7c3e3559e557a239757d41ac02a6 - libraries/Cabal|3dcc425495523ab6142027097cb598a4d2ad810a - libraries/Win32|085b11285b6adbc6484d9c21f5e0b8105556869c - libraries/array|fa295423e7404d3d1d3d82655b2b44d50f045a44 - libraries/base|a57369f54bd25a1de5d477f3c363b3bafd17d168 - libraries/binary|9065df2120254601c68c3a28264fd062abde9354 - libraries/bytestring|caad22630f73e0e9b1b61f4da34f8aefcc6001d8 - libraries/containers|667591b168c804d3eeae503dff1c848ed6852412 - libraries/directory|d44f52926278319286804d8333149dd13d4ecc82 - libraries/dph|b23b45a9e8fce985327b076997d61ab0ddc7b2f7 - libraries/extensible-exceptions|e77722871a5812d52c467e3a8fd9c7b97cdec521 - libraries/filepath|fd381017dca45de5c94dac85a6233516a6b6963d - libraries/ghc-prim|0a84a755e1248b4d50f6535a0ce75af0bb21b0ad - libraries/haskeline|8787a64417500efffc9c48032ee7c37315fb2547 - libraries/haskell2010|470b34b6c0336339eac9fbcfb6020e46b5154bfe - libraries/haskell98|5590c0f042d6d07352e0bf49cedeef5ba0821f23 - libraries/hoopl|b98db91cd0c53ddb2c275c04823f9c379774104b - libraries/hpc|7c726abec939b11af1ecf89740ca8d04e6a1360d - libraries/integer-gmp|65c73842bca2f075b65f418a5ff34753b185e0d7 - libraries/integer-simple|10202111c59f0695ef782d1ec9e6fc992933fc9a - libraries/mtl|a41470c1890073e678f0dca2a9ef4c02d55bf7c6 - libraries/old-locale|104e7e5a7b33424f34e98825a0d4ccb7614ca7c2 - libraries/old-time|81e0c8a4b98d4b084eefe75bedf91a44edd31888 - libraries/pretty|036fb8dfbb9d4a787fcd150c2756b4899be4e942 - libraries/primitive|74479e07b92b8859eae473e5cc86b40decae1d6e - libraries/process|68ba490d6691f55eab19a249379144831055e2ac - libraries/random|3fb0e9e42b54d7b01b794fc27d4d678d7d74ff0e - libraries/template-haskell|02362d12e5ae0af20d637eec97db51f6827a1625 - libraries/terminfo|baec6aff59d13ba294b370f9563e8068706392ce - libraries/unix|f55638fb5c6badd385c51a41de7ff96ef106de42 - libraries/utf8-string|ec2b85940a256dbc8771e5e2065ca8f76cc802d0 - libraries/vector|1e72d46bdc8488a84558b64ac63632cef1d8a695 - libraries/xhtml|cb2cf6c34d815fdf4ed74efeb65e1993e7bda514 - testsuite|26c608a0c31d56917099e4f48bf58c1d1e92e61c - utils/haddock|d54959189f33105ed09a59efee5ba34f53369282 - utils/hsc2hs|f8cbf37ab28ab4512d932678c08c263aa412e008 - - - -First validating in the context of a slightly stale GHC head -(7.3.20110727) on a mac. - - -[2011.09.30] Redoing timings after bugfix in version 1.0.1.1 ------------------------------------------------------------- - -It looks like there has been serious performance regression (3.33ghz -nehalem still). - - How many random numbers can we generate in a second on one thread? - Cost of rdtsc (ffi call): 38 - Approx getCPUTime calls per second: 7,121 - Approx clock frequency: 96,610,524 - First, timing System.Random.next: - 148,133,038 randoms generated [constant zero gen] ~ 0.65 cycles/int - 12,656,455 randoms generated [System.Random stdGen/next] ~ 7.63 cycles/int - - Second, timing System.Random.random at different types: - 676,066 randoms generated [System.Random Ints] ~ 143 cycles/int - 3,917,247 randoms generated [System.Random Word16] ~ 24.66 cycles/int - 2,231,460 randoms generated [System.Random Floats] ~ 43.29 cycles/int - 2,269,993 randoms generated [System.Random CFloats] ~ 42.56 cycles/int - 686,363 randoms generated [System.Random Doubles] ~ 141 cycles/int - 2,165,679 randoms generated [System.Random CDoubles] ~ 44.61 cycles/int - 713,702 randoms generated [System.Random Integers] ~ 135 cycles/int - 3,647,551 randoms generated [System.Random Bools] ~ 26.49 cycles/int - 4,296,919 randoms generated [System.Random Chars] ~ 22.48 cycles/int - - Next timing range-restricted System.Random.randomR: - 4,307,214 randoms generated [System.Random Ints] ~ 22.43 cycles/int - 4,068,982 randoms generated [System.Random Word16s] ~ 23.74 cycles/int - 2,059,264 randoms generated [System.Random Floats] ~ 46.92 cycles/int - 1,960,359 randoms generated [System.Random CFloats] ~ 49.28 cycles/int - 678,978 randoms generated [System.Random Doubles] ~ 142 cycles/int - 2,009,665 randoms generated [System.Random CDoubles] ~ 48.07 cycles/int - 4,296,452 randoms generated [System.Random Integers] ~ 22.49 cycles/int - 3,689,999 randoms generated [System.Random Bools] ~ 26.18 cycles/int - 4,367,577 randoms generated [System.Random Chars] ~ 22.12 cycles/int - 6,650 randoms generated [System.Random BIG Integers] ~ 14,528 cycles/int - diff --git a/LICENSE b/LICENSE index 06bb64148..0c53b4a83 100644 --- a/LICENSE +++ b/LICENSE @@ -1,63 +1,26 @@ -This library (libraries/base) is derived from code from two -sources: - - * Code from the GHC project which is largely (c) The University of - Glasgow, and distributable under a BSD-style license (see below), - - * Code from the Haskell 98 Report which is (c) Simon Peyton Jones - and freely redistributable (but see the full license for - restrictions). - -The full text of these licenses is reproduced below. Both of the -licenses are BSD-style or compatible. - ------------------------------------------------------------------------------ - -The Glasgow Haskell Compiler License - -Copyright 2004, The University Court of the University of Glasgow. +Copyright (c) 2017, Carter Tazio Schonwald 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 name of the University nor the names of its contributors may be -used to endorse or promote products derived from this software without -specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY COURT OF THE UNIVERSITY OF -GLASGOW AND THE 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 -UNIVERSITY COURT OF THE UNIVERSITY OF GLASGOW OR THE 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. - ------------------------------------------------------------------------------ - -Code derived from the document "Report on the Programming Language -Haskell 98", is distributed under the following license: - - Copyright (c) 2002 Simon Peyton Jones - - The authors intend this Report to belong to the entire Haskell - community, and so we grant permission to copy and distribute it for - any purpose, provided that it is reproduced in its entirety, - including this Notice. Modified versions of this Report may also be - copied and distributed for any purpose, provided that the modified - version is clearly presented as such, and that it does not claim to - be a definition of the Haskell 98 Language. - ------------------------------------------------------------------------------ +modification, are permitted provided that the following conditions are +met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +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 deleted file mode 100644 index 47b75896b..000000000 --- a/README.md +++ /dev/null @@ -1,18 +0,0 @@ -The Haskell Standard Library -- Random Number Generation -======================================================== -[![Build Status](https://secure.travis-ci.org/haskell/random.svg?branch=master)](http://travis-ci.org/haskell/random) - -This library provides a basic interface for (splittable) random number generators. - -The API documentation can be found here: - - http://hackage.haskell.org/package/random/docs/System-Random.html - -A module supplying this interface is required for Haskell 98 (but not Haskell -2010). An older [version] -(https://downloads.haskell.org/~ghc/latest/docs/html/libraries/haskell98-2.0.0.3/Random.html) -of this library is included with GHC in the haskell98 package. This newer -version, with compatible api, is included in the [Haskell Platform] -(http://www.haskell.org/platform/contents.html). - -Please report bugs in the Github [issue tracker] (https://github.com/haskell/random/issues) (no longer in the GHC trac). diff --git a/Setup.hs b/Setup.hs index 6fa548caf..9a994af67 100644 --- a/Setup.hs +++ b/Setup.hs @@ -1,6 +1,2 @@ -module Main (main) where - import Distribution.Simple - -main :: IO () main = defaultMain diff --git a/System/Random.hs b/System/Random.hs deleted file mode 100644 index dfd20883b..000000000 --- a/System/Random.hs +++ /dev/null @@ -1,609 +0,0 @@ -#if __GLASGOW_HASKELL__ >= 701 -{-# LANGUAGE Trustworthy #-} -#endif - ------------------------------------------------------------------------------ --- | --- Module : System.Random --- Copyright : (c) The University of Glasgow 2001 --- License : BSD-style (see the file LICENSE in the 'random' repository) --- --- Maintainer : libraries@haskell.org --- Stability : stable --- Portability : portable --- --- This library deals with the common task of pseudo-random number --- generation. The library makes it possible to generate repeatable --- results, by starting with a specified initial random number generator, --- or to get different results on each run by using the system-initialised --- generator or by supplying a seed from some other source. --- --- The library is split into two layers: --- --- * A core /random number generator/ provides a supply of bits. --- The class 'RandomGen' provides a common interface to such generators. --- The library provides one instance of 'RandomGen', the abstract --- data type 'StdGen'. Programmers may, of course, supply their own --- instances of 'RandomGen'. --- --- * The class 'Random' provides a way to extract values of a particular --- type from a random number generator. For example, the 'Float' --- 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. --- ------------------------------------------------------------------------------ - -#include "MachDeps.h" - -module System.Random - ( - - -- $intro - - -- * Random number generators - -#ifdef ENABLE_SPLITTABLEGEN - RandomGen(next, genRange) - , SplittableGen(split) -#else - RandomGen(next, genRange, split) -#endif - -- ** Standard random number generators - , StdGen - , mkStdGen - - -- ** The global random number generator - - -- $globalrng - - , getStdRandom - , getStdGen - , setStdGen - , newStdGen - - -- * Random values of various types - , Random ( random, randomR, - randoms, randomRs, - randomIO, randomRIO ) - - -- * References - -- $references - - ) where - -import Prelude - -import Data.Bits -import Data.Int -import Data.Word -import Foreign.C.Types - -#ifdef __NHC__ -import CPUTime ( getCPUTime ) -import Foreign.Ptr ( Ptr, nullPtr ) -import Foreign.C ( CTime, CUInt ) -#else -import System.CPUTime ( getCPUTime ) -import Data.Time ( getCurrentTime, UTCTime(..) ) -import Data.Ratio ( numerator, denominator ) -#endif -import Data.Char ( isSpace, chr, ord ) -import System.IO.Unsafe ( unsafePerformIO ) -import Data.IORef ( IORef, newIORef, readIORef, writeIORef ) -#if MIN_VERSION_base (4,6,0) -import Data.IORef ( atomicModifyIORef' ) -#else -import Data.IORef ( atomicModifyIORef ) -#endif -import Numeric ( readDec ) - -#ifdef __GLASGOW_HASKELL__ -import GHC.Exts ( build ) -#else --- | A dummy variant of build without fusion. -{-# INLINE build #-} -build :: ((a -> [a] -> [a]) -> [a] -> [a]) -> [a] -build g = g (:) [] -#endif - -#if !MIN_VERSION_base (4,6,0) -atomicModifyIORef' :: IORef a -> (a -> (a,b)) -> IO b -atomicModifyIORef' ref f = do - b <- atomicModifyIORef ref - (\x -> let (a, b) = f x - in (a, a `seq` b)) - b `seq` return b -#endif - --- The standard nhc98 implementation of Time.ClockTime does not match --- the extended one expected in this module, so we lash-up a quick --- replacement here. -#ifdef __NHC__ -foreign import ccall "time.h time" readtime :: Ptr CTime -> IO CTime -getTime :: IO (Integer, Integer) -getTime = do CTime t <- readtime nullPtr; return (toInteger t, 0) -#else -getTime :: IO (Integer, Integer) -getTime = do - utc <- getCurrentTime - let daytime = toRational $ utctDayTime utc - return $ quotRem (numerator daytime) (denominator daytime) -#endif - --- | The class 'RandomGen' provides a common interface to random number --- generators. --- -#ifdef ENABLE_SPLITTABLEGEN --- Minimal complete definition: 'next'. -#else --- Minimal complete definition: 'next' and 'split'. -#endif - -class RandomGen g where - - -- |The 'next' operation returns an 'Int' that is uniformly distributed - -- in the range returned by 'genRange' (including both end points), - -- and a new generator. - next :: g -> (Int, g) - - -- |The 'genRange' operation yields the range of values returned by - -- the generator. - -- - -- It is required that: - -- - -- * If @(a,b) = 'genRange' g@, then @a < b@. - -- - -- * 'genRange' always returns a pair of defined 'Int's. - -- - -- The second condition ensures that 'genRange' cannot examine its - -- argument, and hence the value it returns can be determined only by the - -- instance of 'RandomGen'. That in turn allows an implementation to make - -- a single call to 'genRange' to establish a generator's range, without - -- being concerned that the generator returned by (say) 'next' might have - -- a different range to the generator passed to 'next'. - -- - -- The default definition spans the full range of 'Int'. - genRange :: g -> (Int,Int) - - -- default method - genRange _ = (minBound, maxBound) - -#ifdef ENABLE_SPLITTABLEGEN --- | The class 'SplittableGen' proivides a way to specify a random number --- generator that can be split into two new generators. -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). - 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@. - -In addition, 'reads' may be used to map an arbitrary string (not necessarily one -produced by 'show') onto a value of type 'StdGen'. In general, the 'Read' -instance of 'StdGen' has the following properties: - -* It guarantees to succeed on any string. - -* 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 - -instance RandomGen StdGen where - next = stdNext - genRange _ = stdRange - -#ifdef ENABLE_SPLITTABLEGEN -instance SplittableGen StdGen 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. --} -stdFromString :: String -> (StdGen, String) -stdFromString s = (mkStdGen num, rest) - where (cs, rest) = splitAt 6 s - num = foldl (\a x -> x + 3 * a) 1 (map ord cs) - - -{- | -The function 'mkStdGen' provides an alternative way of producing an initial -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 - -{- -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 - -createStdGen :: Integer -> StdGen -createStdGen s = mkStdGen32 $ fromIntegral s - -{- | -With a source of random number supply in hand, the 'Random' class allows the -programmer to extract random values of a variety of types. - -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. - randomR :: RandomGen g => (a,a) -> g -> (a,g) - - -- | The same as 'randomR', but using a default range determined by the type: - -- - -- * For bounded types (instances of 'Bounded', such as 'Char'), - -- the range is normally the whole type. - -- - -- * For fractional types, the range is normally the semi-closed interval - -- @[0,1)@. - -- - -- * For 'Integer', the range is (arbitrarily) the range of 'Int'. - random :: RandomGen g => g -> (a, g) - - -- | Plural variant of 'randomR', producing an infinite list of - -- random values instead of returning a new generator. - {-# INLINE randomRs #-} - randomRs :: RandomGen g => (a,a) -> g -> [a] - randomRs ival g = build (\cons _nil -> buildRandoms cons (randomR ival) g) - - -- | Plural variant of 'random', producing an infinite list of - -- random values instead of returning a new generator. - {-# INLINE randoms #-} - randoms :: RandomGen g => g -> [a] - randoms g = build (\cons _nil -> buildRandoms cons random g) - - -- | A variant of 'randomR' that uses the global random number generator - -- (see "System.Random#globalrng"). - randomRIO :: (a,a) -> IO a - randomRIO range = getStdRandom (randomR range) - - -- | A variant of 'random' that uses the global random number generator - -- (see "System.Random#globalrng"). - randomIO :: IO a - randomIO = getStdRandom random - --- | Produce an infinite list-equivalent of random values. -{-# INLINE buildRandoms #-} -buildRandoms :: RandomGen g - => (a -> as -> as) -- ^ E.g. '(:)' but subject to fusion - -> (g -> (a,g)) -- ^ E.g. 'random' - -> g -- ^ A 'RandomGen' instance - -> as -buildRandoms cons rand = go - where - -- The seq fixes part of #4218 and also makes fused Core simpler. - go g = x `seq` (x `cons` go g') where (x,g') = rand g - - -instance Random Integer where - randomR ival g = randomIvalInteger ival g - random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g - -instance Random Int where randomR = randomIvalIntegral; random = randomBounded -instance Random Int8 where randomR = randomIvalIntegral; random = randomBounded -instance Random Int16 where randomR = randomIvalIntegral; random = randomBounded -instance Random Int32 where randomR = randomIvalIntegral; random = randomBounded -instance Random Int64 where randomR = randomIvalIntegral; random = randomBounded - -#ifndef __NHC__ --- Word is a type synonym in nhc98. -instance Random Word where randomR = randomIvalIntegral; random = randomBounded -#endif -instance Random Word8 where randomR = randomIvalIntegral; random = randomBounded -instance Random Word16 where randomR = randomIvalIntegral; random = randomBounded -instance Random Word32 where randomR = randomIvalIntegral; random = randomBounded -instance Random Word64 where randomR = randomIvalIntegral; random = randomBounded - -instance Random CChar where randomR = randomIvalIntegral; random = randomBounded -instance Random CSChar where randomR = randomIvalIntegral; random = randomBounded -instance Random CUChar where randomR = randomIvalIntegral; random = randomBounded -instance Random CShort where randomR = randomIvalIntegral; random = randomBounded -instance Random CUShort where randomR = randomIvalIntegral; random = randomBounded -instance Random CInt where randomR = randomIvalIntegral; random = randomBounded -instance Random CUInt where randomR = randomIvalIntegral; random = randomBounded -instance Random CLong where randomR = randomIvalIntegral; random = randomBounded -instance Random CULong where randomR = randomIvalIntegral; random = randomBounded -instance Random CPtrdiff where randomR = randomIvalIntegral; random = randomBounded -instance Random CSize where randomR = randomIvalIntegral; random = randomBounded -instance Random CWchar where randomR = randomIvalIntegral; random = randomBounded -instance Random CSigAtomic where randomR = randomIvalIntegral; random = randomBounded -instance Random CLLong where randomR = randomIvalIntegral; random = randomBounded -instance Random CULLong where randomR = randomIvalIntegral; random = randomBounded -instance Random CIntPtr where randomR = randomIvalIntegral; random = randomBounded -instance Random CUIntPtr where randomR = randomIvalIntegral; random = randomBounded -instance Random CIntMax where randomR = randomIvalIntegral; random = randomBounded -instance Random CUIntMax where randomR = randomIvalIntegral; random = randomBounded - -instance Random Char where - randomR (a,b) g = - case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of - (x,g') -> (chr x, g') - random g = randomR (minBound,maxBound) g - -instance Random Bool where - randomR (a,b) g = - case (randomIvalInteger (bool2Int a, bool2Int b) g) of - (x, g') -> (int2Bool x, g') - where - bool2Int :: Bool -> Integer - bool2Int False = 0 - bool2Int True = 1 - - int2Bool :: Int -> Bool - int2Bool 0 = False - int2Bool _ = True - - random g = randomR (minBound,maxBound) g - -{-# INLINE randomRFloating #-} -randomRFloating :: (Fractional a, Num a, Ord a, Random a, RandomGen g) => (a, a) -> g -> (a, g) -randomRFloating (l,h) g - | l>h = randomRFloating (h,l) g - | otherwise = let (coef,g') = random g in - (2.0 * (0.5*l + coef * (0.5*h - 0.5*l)), g') -- avoid overflow - -instance Random Double where - randomR = randomRFloating - random rng = - case random rng of - (x,rng') -> - -- We use 53 bits of randomness corresponding to the 53 bit significand: - ((fromIntegral (mask53 .&. (x::Int64)) :: Double) - / fromIntegral twoto53, rng') - where - twoto53 = (2::Int64) ^ (53::Int64) - mask53 = twoto53 - 1 - -instance Random Float where - randomR = randomRFloating - random rng = - -- TODO: Faster to just use 'next' IF it generates enough bits of randomness. - case random rng of - (x,rng') -> - -- We use 24 bits of randomness corresponding to the 24 bit significand: - ((fromIntegral (mask24 .&. (x::Int32)) :: Float) - / fromIntegral twoto24, rng') - -- Note, encodeFloat is another option, but I'm not seeing slightly - -- worse performance with the following [2011.06.25]: --- (encodeFloat rand (-24), rng') - where - mask24 = twoto24 - 1 - twoto24 = (2::Int32) ^ (24::Int32) - --- CFloat/CDouble are basically the same as a Float/Double: -instance Random CFloat where - randomR = randomRFloating - random rng = case random rng of - (x,rng') -> (realToFrac (x::Float), rng') - -instance Random CDouble where - randomR = randomRFloating - -- A MYSTERY: - -- Presently, this is showing better performance than the Double instance: - -- (And yet, if the Double instance uses randomFrac then its performance is much worse!) - random = randomFrac - -- 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) - --- The two integer functions below take an [inclusive,inclusive] range. -randomIvalIntegral :: (RandomGen g, Integral a) => (a, a) -> g -> (a, g) -randomIvalIntegral (l,h) = randomIvalInteger (toInteger l, toInteger h) - -{-# SPECIALIZE randomIvalInteger :: (Num a) => - (Integer, Integer) -> StdGen -> (a, StdGen) #-} - -randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g) -randomIvalInteger (l,h) rng - | l > h = randomIvalInteger (h,l) rng - | otherwise = case (f 1 0 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng') - where - (genlo, genhi) = genRange rng - b = fromIntegral genhi - fromIntegral genlo + 1 - - -- Probabilities of the most likely and least likely result - -- will differ at most by a factor of (1 +- 1/q). Assuming the RandomGen - -- is uniform, of course - - -- On average, log q / log b more random values will be generated - -- than the minimum - q = 1000 - k = h - l + 1 - magtgt = k * q - - -- generate random values until we exceed the target magnitude - f mag v g | mag >= magtgt = (v, g) - | otherwise = v' `seq`f (mag*b) v' g' where - (x,g') = next g - v' = (v * b + (fromIntegral x - fromIntegral genlo)) - - --- The continuous functions on the other hand take an [inclusive,exclusive) range. -randomFrac :: (RandomGen g, Fractional a) => g -> (a, g) -randomFrac = randomIvalDouble (0::Double,1) realToFrac - -randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g) -randomIvalDouble (l,h) fromDouble rng - | l > h = randomIvalDouble (h,l) fromDouble rng - | otherwise = - case (randomIvalInteger (toInteger (minBound::Int32), toInteger (maxBound::Int32)) rng) of - (x, rng') -> - let - scaled_x = - fromDouble (0.5*l + 0.5*h) + -- previously (l+h)/2, overflowed - fromDouble ((0.5*h - 0.5*l) / (0.5 * realToFrac int32Count)) * -- avoid overflow - fromIntegral (x::Int32) - in - (scaled_x, 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# - -There is a single, implicit, global random number generator of type -'StdGen', held in some global variable maintained by the 'IO' monad. It is -initialised automatically in some system-dependent fashion, for example, by -using the time of day, or Linux's kernel random number generator. To get -deterministic behaviour, use 'setStdGen'. --} - --- |Sets the global random number generator. -setStdGen :: StdGen -> IO () -setStdGen sgen = writeIORef theStdGen sgen - --- |Gets the global random number generator. -getStdGen :: IO StdGen -getStdGen = readIORef theStdGen - -theStdGen :: IORef StdGen -theStdGen = unsafePerformIO $ do - rng <- mkStdRNG 0 - newIORef rng - --- |Applies 'split' to the current global random generator, --- updates it with one of the results, and returns the other. -newStdGen :: IO StdGen -newStdGen = atomicModifyIORef' theStdGen split - -{- |Uses the supplied function to get a value from the current global -random generator, and updates the global generator with the new generator -returned by the function. For example, @rollDice@ gets a random integer -between 1 and 6: - -> rollDice :: IO Int -> rollDice = getStdRandom (randomR (1,6)) - --} - -getStdRandom :: (StdGen -> (a,StdGen)) -> IO a -getStdRandom f = atomicModifyIORef' theStdGen (swap . f) - where swap (v,g) = (g,v) - -{- $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, -, 1998. - -5. Pierre #LEcuyer# L'Ecuyer, /Efficient and portable combined random -number generators/, Comm ACM, 31(6), Jun 1988, pp742-749. - -The Web site is a great source of information. - --} diff --git a/prologue.txt b/prologue.txt deleted file mode 100644 index ea0344b9c..000000000 --- a/prologue.txt +++ /dev/null @@ -1 +0,0 @@ -Random number library. diff --git a/random.cabal b/random.cabal index a28063c90..ac22c6816 100644 --- a/random.cabal +++ b/random.cabal @@ -1,70 +1,70 @@ -name: random -version: 1.1 - - - - -license: BSD3 -license-file: LICENSE -maintainer: core-libraries-committee@haskell.org -bug-reports: https://github.com/haskell/random/issues -synopsis: random number library -category: System -description: - This package provides a basic random number generation - library, including the ability to split random number - generators. - -extra-source-files: - .travis.yml - README.md - CHANGELOG.md - .gitignore - .darcs-boring - - - -build-type: Simple --- cabal-version 1.8 needed because "the field 'build-depends: random' refers --- to a library which is defined within the same package" -cabal-version: >= 1.8 - - - -Library - exposed-modules: - System.Random - extensions: CPP - GHC-Options: -O2 - build-depends: base >= 3 && < 5, time - -source-repository head - type: git - location: https://github.com/haskell/random.git - --- To run the Test-Suite: --- $ cabal configure --enable-tests --- $ cabal test --show-details=always --test-options="+RTS -M1M -RTS" - -Test-Suite T7936 - type: exitcode-stdio-1.0 - main-is: T7936.hs - hs-source-dirs: tests - build-depends: base >= 3 && < 5, random - ghc-options: -rtsopts -O2 - -Test-Suite TestRandomRs - type: exitcode-stdio-1.0 - main-is: TestRandomRs.hs - hs-source-dirs: tests - build-depends: base >= 3 && < 5, random - ghc-options: -rtsopts -O2 - -- TODO. Why does the following not work? - --test-options: +RTS -M1M -RTS - -Test-Suite TestRandomIOs - type: exitcode-stdio-1.0 - main-is: TestRandomIOs.hs - hs-source-dirs: tests - build-depends: base >= 3 && < 5, random - ghc-options: -rtsopts -O2 +-- 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) +-- for standards guiding when and how versions should be incremented. +-- https://wiki.haskell.org/Package_versioning_policy +-- PVP summary: +-+------- breaking API changes +-- | | +----- non-breaking API additions +-- | | | +--- code changes with no API change +version: 2.0.0.0 + +-- A short (one-line) description of the package. +synopsis: Random number generation library for haskell + +-- A longer description of the package. +-- description: + +-- URL for the project homepage or repository. +homepage: http://github.com/cartazio/random + +-- The license under which the package is released. +license: BSD2 + +-- The file containing the license text. +license-file: LICENSE + +-- The package author(s). +author: Carter Tazio Schonwald + +-- An email address to which users can send suggestions, bug reports, and +-- patches. +maintainer: carter at wellposed dot com + +-- A copyright notice. +-- copyright: + +category: Math + +build-type: Simple + +-- Extra files to be distributed with the package, such as examples or a +-- README. +extra-source-files: ChangeLog.md + +-- Constraint on the version of Cabal needed to build this package. +cabal-version: >=1.10 + + +library + -- Modules exported by the library. + -- exposed-modules: + + -- Modules included in this library but not exported. + -- other-modules: + + -- LANGUAGE extensions used by modules in this package. + -- other-extensions: + + -- Other library packages from which modules are imported. + build-depends: base >=4.9 && <4.10 + + -- Directories containing source files. + hs-source-dirs: src + + -- Base language which the package is written in. + default-language: Haskell2010 + diff --git a/tests/Makefile b/tests/Makefile deleted file mode 100644 index 39c71491b..000000000 --- a/tests/Makefile +++ /dev/null @@ -1,14 +0,0 @@ -# This Makefile runs the tests using GHC's testsuite framework. It -# assumes the package is part of a GHC build tree with the testsuite -# installed in ../../../testsuite. - -TOP=../../../testsuite -include $(TOP)/mk/boilerplate.mk -include $(TOP)/mk/test.mk - - -# Build tests locally without the central GHC testing infrastructure: -local: - ghc --make rangeTest.hs - ghc --make random1283.hs - diff --git a/tests/T7936.hs b/tests/T7936.hs deleted file mode 100644 index cfea911bc..000000000 --- a/tests/T7936.hs +++ /dev/null @@ -1,14 +0,0 @@ --- Test for ticket #7936: --- https://ghc.haskell.org/trac/ghc/ticket/7936 --- --- Used to fail with: --- --- $ cabal test T7936 --test-options="+RTS -M1M -RTS" --- T7936: Heap exhausted; - -module Main where - -import System.Random (newStdGen) -import Control.Monad (replicateM_) - -main = replicateM_ 100000 newStdGen diff --git a/tests/TestRandomIOs.hs b/tests/TestRandomIOs.hs deleted file mode 100644 index d8a00cc7c..000000000 --- a/tests/TestRandomIOs.hs +++ /dev/null @@ -1,20 +0,0 @@ --- Test for ticket #4218 (TestRandomIOs): --- https://ghc.haskell.org/trac/ghc/ticket/4218 --- --- Used to fail with: --- --- $ cabal test TestRandomIOs --test-options="+RTS -M1M -RTS" --- TestRandomIOs: Heap exhausted; - -module Main where - -import Control.Monad (replicateM) -import System.Random (randomIO) - --- Build a list of 5000 random ints in memory (IO Monad is strict), and print --- the last one. --- Should use less than 1Mb of heap space, or we are generating a list of --- unevaluated thunks. -main = do - rs <- replicateM 5000 randomIO :: IO [Int] - print $ last rs diff --git a/tests/TestRandomRs.hs b/tests/TestRandomRs.hs deleted file mode 100644 index cdae106a6..000000000 --- a/tests/TestRandomRs.hs +++ /dev/null @@ -1,22 +0,0 @@ --- Test for ticket #4218 (TestRandomRs): --- https://ghc.haskell.org/trac/ghc/ticket/4218 --- --- Fixed together with ticket #8704 --- https://ghc.haskell.org/trac/ghc/ticket/8704 --- Commit 4695ffa366f659940369f05e419a4f2249c3a776 --- --- Used to fail with: --- --- $ cabal test TestRandomRs --test-options="+RTS -M1M -RTS" --- TestRandomRs: Heap exhausted; - -module Main where - -import Control.Monad (liftM, replicateM) -import System.Random (randomRs, getStdGen) - --- Return the five-thousandth random number: --- Should run in constant space (< 1Mb heap). -main = do - n <- (last . take 5000 . randomRs (0, 1000000)) `liftM` getStdGen - print (n::Integer) diff --git a/tests/all.T b/tests/all.T deleted file mode 100644 index f85ec767a..000000000 --- a/tests/all.T +++ /dev/null @@ -1,11 +0,0 @@ -setTestOpts(reqlib('random')) - -test('rangeTest', - normal, - compile_and_run, - ['']) - -test('random1283', - reqlib('containers'), - compile_and_run, - [' -package containers']) diff --git a/tests/random1283.hs b/tests/random1283.hs deleted file mode 100644 index 33fc48862..000000000 --- a/tests/random1283.hs +++ /dev/null @@ -1,40 +0,0 @@ -import Control.Concurrent -import Control.Monad hiding (empty) -import Data.Sequence (ViewL(..), empty, fromList, viewl, (<|), (|>), (><)) -import System.Random - --- This test - -threads = 4 -samples = 5000 - -main = loopTest threads samples - -loopTest t s = do - isClean <- testRace t s - when (not isClean) $ putStrLn "race condition!" - -testRace t s = do - ref <- liftM (take (t*s) . randoms) getStdGen - iss <- threadRandoms t s - return (isInterleavingOf (ref::[Int]) iss) - -threadRandoms :: Random a => Int -> Int -> IO [[a]] -threadRandoms t s = do - vs <- sequence $ replicate t $ do - v <- newEmptyMVar - forkIO (sequence (replicate s randomIO) >>= putMVar v) - return v - mapM takeMVar vs - -isInterleavingOf xs yss = iio xs (viewl $ fromList yss) EmptyL where - iio (x:xs) ((y:ys) :< yss) zss - | x /= y = iio (x:xs) (viewl yss) (viewl (fromViewL zss |> (y:ys))) - | x == y = iio xs (viewl ((ys <| yss) >< fromViewL zss)) EmptyL - iio xs ([] :< yss) zss = iio xs (viewl yss) zss - iio [] EmptyL EmptyL = True - iio _ _ _ = False - -fromViewL (EmptyL) = empty -fromViewL (x :< xs) = x <| xs - diff --git a/tests/rangeTest.hs b/tests/rangeTest.hs deleted file mode 100644 index ac62c71dd..000000000 --- a/tests/rangeTest.hs +++ /dev/null @@ -1,135 +0,0 @@ -{-# LANGUAGE CPP #-} -import Control.Monad -import System.Random -import Data.Int -import Data.Word -import Data.Bits -import Foreign.C.Types - --- Take many measurements and record the max/min/average random values. -approxBounds :: (RandomGen g, Random a, Ord a, Num a) => - (g -> (a,g)) -> Int -> a -> (a,a) -> g -> ((a,a,a),g) --- Here we do a little hack to essentiall pass in the type in the last argument: -approxBounds nxt iters unused (explo,exphi) initrng = - if False - then ((unused,unused,unused),undefined) --- else loop initrng iters 100 (-100) 0 -- Oops, can't use minBound/maxBound here. - else loop initrng iters exphi explo 0 -- Oops, can't use minBound/maxBound here. - where - loop rng 0 mn mx sum = ((mn,mx,sum),rng) - loop rng n mn mx sum = - case nxt rng of - (x, rng') -> loop rng' (n-1) (min x mn) (max x mx) (x+sum) - - --- We check that: --- (1) all generated numbers are in bounds --- (2) we get "close" to the bounds --- The with (2) is that we do enough trials to ensure that we can at --- least hit the 90% mark. -checkBounds:: (Real a, Show a, Ord a) => - String -> (Bool, a, a) -> ((a,a) -> StdGen -> ((a, a, t), StdGen)) -> IO () -checkBounds msg (exclusive,lo,hi) fun = - -- (lo,hi) is [inclusive,exclusive) - do putStr$ msg --- ++ ", expected range " ++ show (lo,hi) - ++ ": " - (mn,mx,sum) <- getStdRandom (fun (lo,hi)) - when (mn < lo)$ error$ "broke lower bound: " ++ show mn - when (mx > hi) $ error$ "broke upper bound: " ++ show mx - when (exclusive && mx >= hi)$ error$ "hit upper bound: " ++ show mx - - let epsilon = 0.1 * (toRational hi - toRational lo) - - when (toRational (hi - mx) > epsilon)$ error$ "didn't get close enough to upper bound: "++ show mx - when (toRational (mn - lo) > epsilon)$ error$ "didn't get close enough to lower bound: "++ show mn - putStrLn "Passed" - -boundedRange :: (Num a, Bounded a) => (Bool, a, a) -boundedRange = ( False, minBound, maxBound ) - -trials = 5000 - --- Keep in mind here that on some architectures (e.g. ARM) CChar, CWchar, and CSigAtomic --- are unsigned - -main = - do - checkBounds "Int" boundedRange (approxBounds random trials (undefined::Int)) - checkBounds "Integer" (False, fromIntegral (minBound::Int), fromIntegral (maxBound::Int)) - (approxBounds random trials (undefined::Integer)) - checkBounds "Int8" boundedRange (approxBounds random trials (undefined::Int8)) - checkBounds "Int16" boundedRange (approxBounds random trials (undefined::Int16)) - checkBounds "Int32" boundedRange (approxBounds random trials (undefined::Int32)) - checkBounds "Int64" boundedRange (approxBounds random trials (undefined::Int64)) - checkBounds "Word" boundedRange (approxBounds random trials (undefined::Word)) - checkBounds "Word8" boundedRange (approxBounds random trials (undefined::Word8)) - checkBounds "Word16" boundedRange (approxBounds random trials (undefined::Word16)) - checkBounds "Word32" boundedRange (approxBounds random trials (undefined::Word32)) - checkBounds "Word64" boundedRange (approxBounds random trials (undefined::Word64)) - checkBounds "Double" (True,0.0,1.0) (approxBounds random trials (undefined::Double)) - checkBounds "Float" (True,0.0,1.0) (approxBounds random trials (undefined::Float)) - - checkBounds "CChar" boundedRange (approxBounds random trials (undefined:: CChar)) - checkBounds "CSChar" boundedRange (approxBounds random trials (undefined:: CSChar)) - checkBounds "CUChar" boundedRange (approxBounds random trials (undefined:: CUChar)) - checkBounds "CShort" boundedRange (approxBounds random trials (undefined:: CShort)) - checkBounds "CUShort" boundedRange (approxBounds random trials (undefined:: CUShort)) - checkBounds "CInt" boundedRange (approxBounds random trials (undefined:: CInt)) - checkBounds "CUInt" boundedRange (approxBounds random trials (undefined:: CUInt)) - checkBounds "CLong" boundedRange (approxBounds random trials (undefined:: CLong)) - checkBounds "CULong" boundedRange (approxBounds random trials (undefined:: CULong)) - checkBounds "CPtrdiff" boundedRange (approxBounds random trials (undefined:: CPtrdiff)) - checkBounds "CSize" boundedRange (approxBounds random trials (undefined:: CSize)) - checkBounds "CWchar" boundedRange (approxBounds random trials (undefined:: CWchar)) - checkBounds "CSigAtomic" boundedRange (approxBounds random trials (undefined:: CSigAtomic)) - checkBounds "CLLong" boundedRange (approxBounds random trials (undefined:: CLLong)) - checkBounds "CULLong" boundedRange (approxBounds random trials (undefined:: CULLong)) - checkBounds "CIntPtr" boundedRange (approxBounds random trials (undefined:: CIntPtr)) - checkBounds "CUIntPtr" boundedRange (approxBounds random trials (undefined:: CUIntPtr)) - checkBounds "CIntMax" boundedRange (approxBounds random trials (undefined:: CIntMax)) - checkBounds "CUIntMax" boundedRange (approxBounds random trials (undefined:: CUIntMax)) - - -- Then check all the range-restricted versions: - checkBounds "Int R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined::Int)) - checkBounds "Integer R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined::Integer)) - checkBounds "Int8 R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined::Int8)) - checkBounds "Int8 Rsmall" (False,-50,50) (approxBounds (randomR (-50,50)) trials (undefined::Int8)) - checkBounds "Int8 Rmini" (False,3,4) (approxBounds (randomR (3,4)) trials (undefined::Int8)) - checkBounds "Int8 Rtrivial" (False,3,3) (approxBounds (randomR (3,3)) trials (undefined::Int8)) - - checkBounds "Int16 R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined::Int16)) - checkBounds "Int32 R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined::Int32)) - checkBounds "Int64 R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined::Int64)) - checkBounds "Word R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined::Word)) - checkBounds "Word8 R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined::Word8)) - checkBounds "Word16 R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined::Word16)) - checkBounds "Word32 R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined::Word32)) - checkBounds "Word64 R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined::Word64)) - checkBounds "Double R" (True,10.0,77.0) (approxBounds (randomR (10,77)) trials (undefined::Double)) - checkBounds "Float R" (True,10.0,77.0) (approxBounds (randomR (10,77)) trials (undefined::Float)) - - checkBounds "CChar R" (False,0,100) (approxBounds (randomR (0,100)) trials (undefined:: CChar)) - checkBounds "CSChar R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CSChar)) - checkBounds "CUChar R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CUChar)) - checkBounds "CShort R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CShort)) - checkBounds "CUShort R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CUShort)) - checkBounds "CInt R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CInt)) - checkBounds "CUInt R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CUInt)) - checkBounds "CLong R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CLong)) - checkBounds "CULong R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CULong)) - checkBounds "CPtrdiff R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CPtrdiff)) - checkBounds "CSize R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CSize)) - checkBounds "CWchar R" (False,0,100) (approxBounds (randomR (0,100)) trials (undefined:: CWchar)) - checkBounds "CSigAtomic R" (False,0,100) (approxBounds (randomR (0,100)) trials (undefined:: CSigAtomic)) - checkBounds "CLLong R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CLLong)) - checkBounds "CULLong R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CULLong)) - checkBounds "CIntPtr R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CIntPtr)) - checkBounds "CUIntPtr R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CUIntPtr)) - checkBounds "CIntMax R" (False,-100,100) (approxBounds (randomR (-100,100)) trials (undefined:: CIntMax)) - checkBounds "CUIntMax R" (False,0,200) (approxBounds (randomR (0,200)) trials (undefined:: CUIntMax)) - --- Untested: --- instance Random Char where --- instance Random Bool where --- instance Random Integer where diff --git a/tests/rangeTest.stdout b/tests/rangeTest.stdout deleted file mode 100644 index 55ccaffb4..000000000 --- a/tests/rangeTest.stdout +++ /dev/null @@ -1,67 +0,0 @@ -Int: Passed -Integer: Passed -Int8: Passed -Int16: Passed -Int32: Passed -Int64: Passed -Word: Passed -Word8: Passed -Word16: Passed -Word32: Passed -Word64: Passed -Double: Passed -Float: Passed -CChar: Passed -CSChar: Passed -CUChar: Passed -CShort: Passed -CUShort: Passed -CInt: Passed -CUInt: Passed -CLong: Passed -CULong: Passed -CPtrdiff: Passed -CSize: Passed -CWchar: Passed -CSigAtomic: Passed -CLLong: Passed -CULLong: Passed -CIntPtr: Passed -CUIntPtr: Passed -CIntMax: Passed -CUIntMax: Passed -Int R: Passed -Integer R: Passed -Int8 R: Passed -Int8 Rsmall: Passed -Int8 Rmini: Passed -Int8 Rtrivial: Passed -Int16 R: Passed -Int32 R: Passed -Int64 R: Passed -Word R: Passed -Word8 R: Passed -Word16 R: Passed -Word32 R: Passed -Word64 R: Passed -Double R: Passed -Float R: Passed -CChar R: Passed -CSChar R: Passed -CUChar R: Passed -CShort R: Passed -CUShort R: Passed -CInt R: Passed -CUInt R: Passed -CLong R: Passed -CULong R: Passed -CPtrdiff R: Passed -CSize R: Passed -CWchar R: Passed -CSigAtomic R: Passed -CLLong R: Passed -CULLong R: Passed -CIntPtr R: Passed -CUIntPtr R: Passed -CIntMax R: Passed -CUIntMax R: Passed diff --git a/tests/slowness.hs b/tests/slowness.hs deleted file mode 100644 index 162a4cff0..000000000 --- a/tests/slowness.hs +++ /dev/null @@ -1,55 +0,0 @@ - --- http://hackage.haskell.org/trac/ghc/ticket/427 - --- Two (performance) problems in one: - -{-# OPTIONS -fffi #-} -module Main (main) where - -import Control.Monad -import System.Random - -foreign import ccall unsafe "random" _crandom :: IO Int -foreign import ccall unsafe "stdlib.hs" rand :: IO Int - -randomInt :: (Int, Int) -> IO Int -randomInt (min,max) = do --- n <- _crandom - n <- rand - return $ min + n `rem` range - where - range = max - min + 1 - -main = replicateM_ (5*10^6) $ do - x <- randomRIO (0::Int,1000) :: IO Int --- x <- randomInt (0::Int,1000) :: IO Int - x `seq` return () - return () - --- First, without the "seq" at the end, hardly anything is --- evaluated and we're building huge amounts of thunks. --- Three ideas about this one: --- - Blame the user :) --- - data StdGen = StdGen !Int !Int --- Use strict fields in StdGen. Doesn't actually help --- (at least in this example). --- - Force evaluation of the StdGen in getStdRandom. --- Does help in this example, but also changes behaviour --- of the library: --- x <- randomRIO undefined --- currently dies only when x (or the result of a later --- randomRIO) is evaluated. This change causes it to die --- immediately. - --- Second, even _with_ the "seq", replacing "randomRIO" by --- "randomInt" speeds the thing up with a factor of about --- 30. (2 to 3.6, in a "real world" university practicum --- exercise of 900 lines of code) --- Even given the fact that they're not really doing the --- same thing, this seems rather much :( - --------------------------------------------------------------------------------- - --- [2011.06.28] RRN: --- I'm currently seeing 1425 ms vs 43 ms for the above. 33X --- difference. If I use rand() instead it's about 52ms. From 3e6efad4f6e069e62cbf4e7eafc5d62ef9848eb0 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sun, 15 Jan 2017 10:30:12 -0500 Subject: [PATCH 08/31] pure Splitmix implementation --- random.cabal | 34 +++--- src/System/Random/SplitMix/Internal.hs | 156 +++++++++++++++++++++++++ 2 files changed, 174 insertions(+), 16 deletions(-) create mode 100644 src/System/Random/SplitMix/Internal.hs diff --git a/random.cabal b/random.cabal index ac22c6816..6f023fdee 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 000000000..e585586a8 --- /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; +-} From c0fb9111e1a8c6ada51eb8e9630b74ec8c7d2bf0 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Fri, 10 Feb 2017 14:15:22 -0500 Subject: [PATCH 09/31] initial attempt to see how far i could get with s -> (# m a , s #) monad flavor thats impossible, so switchign to s -> m (a,s) flavor --- random.cabal | 2 +- src/System/Random/SplitMix/Internal.hs | 129 ++++++++++++------------- 2 files changed, 62 insertions(+), 69 deletions(-) diff --git a/random.cabal b/random.cabal index 6f023fdee..2bc6e1fa6 100644 --- a/random.cabal +++ b/random.cabal @@ -60,7 +60,7 @@ library -- other-extensions: -- Other library packages from which modules are imported. - build-depends: base >=4.8 && <4.10 + build-depends: base >=4.8 && <4.10, ghc-prim -- Directories containing source files. hs-source-dirs: src diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index e585586a8..198393eb5 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -1,8 +1,14 @@ -{-# LANGUAGE ScopedTypeVariables, BangPatterns, UnboxedTuples #-} +{-# LANGUAGE ScopedTypeVariables, BangPatterns, UnboxedTuples, MagicHash, GADTs #-} +{-# LANGUAGE DeriveFunctor #-} + module System.Random.SplitMix.Internal( - --mix32, - xorShiftR + nextSeedSplitMix + ,splitGeneratorSplitMix + ,nextWord64SplitMix + ,SplitMix64(..) + ,Random(..) + ,RandomT(..) ) where import qualified Data.Bits as DB @@ -56,12 +62,19 @@ mixGamma = \ w -> runIdentity $! then return (mixedGamma `xor` 0xaaaaaaaaaaaaaaaa) else return mixedGamma +{- +theres a few different alternatives we could do for the RNG state + +-- this isn't quite expressible +type SplitMix64 = (# Word64# , Word64# #) +-} data SplitMix64 = SplitMix64 { sm64seed :: {-# UNPACK #-} !Word64 ,sm64Gamma :: {-# UNPACK #-} !Word64 } + advanceSplitMix :: SplitMix64 -> SplitMix64 advanceSplitMix (SplitMix64 sd gamma) = SplitMix64 (sd + gamma) gamma @@ -71,86 +84,66 @@ nextSeedSplitMix gen@(SplitMix64 result _) = newgen `seq` (# result,newgen #) newgen = advanceSplitMix gen -newtype RandomM a = RandomM (SplitMix64 -> (# a , SplitMix64 #)) +newtype Random a = Random# (SplitMix64 -> (# a , SplitMix64 #)) + --deriving Functor -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; -} +newtype RandomT m a = RandomT# { unRandomT# :: (SplitMix64 -> (# m a , SplitMix64 #)) } -inline void advance(struct SplitMix64* generator); -inline uint64_t next_seed(struct SplitMix64* generator); +instance Functor m => Functor (RandomT m) where + fmap = \ f (RandomT# mf) -> + RandomT# $ \ seed -> + let (# !ma , !s' #) = mf seed + !mb = fmap f ma + in (# mb , s' #) -inline void advance(struct SplitMix64* generator) { - generator->seed += generator->gamma; -} +instance Applicative m => Applicative (RandomT m) where + pure = \ x -> RandomT# $ \ s -> (# pure x , s #) + (<*>) = \ (RandomT# frmb) (RandomT# rma) -> RandomT# $ \ s -> + let (# !fseed, !maseed #) = splitGeneratorSplitMix s + (# !mf , _boringSeed #) = frmb fseed + (# !ma , newSeed #) = rma maseed + in (# mf <*> ma , newSeed #) -inline uint64_t next_seed(struct SplitMix64* generator) { - uint64_t result = generator->seed; - advance(generator); - return result; -} +instance Monad m => Monad (RandomT m) where + (>>=) = \ (RandomT# ma) mf -> + RandomT# $ \ s -> + let + (# splitSeed, nextSeed #) = splitGeneratorSplitMix s + (# maRes, _boringSeed #) = ma splitSeed + (# mfRes , resultSeed #) -uint64_t next_int64(struct SplitMix64* generator) { - return mix64(next_seed(generator)); -} +{- +there are two models of RandomT m a we could do -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; - } - } -} +1) s -> (m a , s) +or +2) s -> m (a,s) -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; -} +-- The 'return' function leaves the state unchanged, while @>>=@ uses +-- split on the rng state so that the final state of the first computation +-- is independent of the second ... +so lets try writing an instance using 1 +-} -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))); +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 -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; --} From 644c184b95c2a58734116c83350364ca7d85b81e Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Fri, 10 Feb 2017 14:46:57 -0500 Subject: [PATCH 10/31] RandomT in the style of s -> m (a,s) is working --- src/System/Random/SplitMix/Internal.hs | 63 +++++++++++++++++++------- 1 file changed, 46 insertions(+), 17 deletions(-) diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index 198393eb5..6c6dac872 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -84,36 +84,65 @@ nextSeedSplitMix gen@(SplitMix64 result _) = newgen `seq` (# result,newgen #) newgen = advanceSplitMix gen -newtype Random a = Random# (SplitMix64 -> (# a , SplitMix64 #)) - --deriving Functor +newtype Random a = Random# {unRandom# :: SplitMix64 -> (# a , SplitMix64 #)} +instance Functor Random where + fmap = \ f (Random# mf) -> + Random# $ \ seed -> + let (# !a , !s' #) = mf seed + !b = f a + in (# b , s' #) +instance Applicative Random where + pure = \ x -> Random# $ \ s -> (# x , s #) + (<*>) = \ (Random# frmb) (Random# rma) -> Random# $ \ s -> + let (# fseed, maseed #) = splitGeneratorSplitMix s + (# f , _boringSeed #) = frmb fseed + (# a , newSeed #) = rma maseed + in (# f a , newSeed #) -newtype RandomT m a = RandomT# { unRandomT# :: (SplitMix64 -> (# m a , SplitMix64 #)) } +instance Monad Random where + (>>=) = + \(Random# ma) f -> + Random# $ \ s -> + let (# splitSeed , nextSeed #) = splitGeneratorSplitMix s + (# a, _boringSeed #) = ma splitSeed + in unRandom# (f a) nextSeed + + + +newtype RandomT m a = RandomT# { unRandomT# :: (SplitMix64 -> m (a , SplitMix64) ) } instance Functor m => Functor (RandomT m) where fmap = \ f (RandomT# mf) -> RandomT# $ \ seed -> - let (# !ma , !s' #) = mf seed - !mb = fmap f ma - in (# mb , s' #) + fmap (\(a,s) -> (f a, s) ) $ mf seed instance Applicative m => Applicative (RandomT m) where - pure = \ x -> RandomT# $ \ s -> (# pure x , s #) + pure = \ x -> RandomT# $ \ s -> pure ( x , s ) (<*>) = \ (RandomT# frmb) (RandomT# rma) -> RandomT# $ \ s -> let (# !fseed, !maseed #) = splitGeneratorSplitMix s - (# !mf , _boringSeed #) = frmb fseed - (# !ma , newSeed #) = rma maseed - in (# mf <*> ma , newSeed #) - + mfOldSeed= frmb fseed + mArgNewSeed = rma maseed + in (fmap (\(f,_s)-> \(x,newSeed)-> (f x, newSeed) ) mfOldSeed) + <*> mArgNewSeed instance Monad m => Monad (RandomT m) where - (>>=) = \ (RandomT# ma) mf -> - RandomT# $ \ s -> - let - (# splitSeed, nextSeed #) = splitGeneratorSplitMix s - (# maRes, _boringSeed #) = ma splitSeed - (# mfRes , resultSeed #) + + (>>=) = \ (RandomT# ma) f -> RandomT# $ \ s -> + let (# fseed, nextSeed #) = splitGeneratorSplitMix s + in + do + (a,_boring) <- ma fseed + unRandomT# (f a) nextSeed + +--instance Monad m => Monad (RandomT m) where +-- (>>=) = \ (RandomT# ma) mf -> +-- RandomT# $ \ s -> +-- let +-- (# splitSeed, nextSeed #) = splitGeneratorSplitMix s +-- (# maRes, _boringSeed #) = ma splitSeed +-- (# mfRes , resultSeed #) {- there are two models of RandomT m a we could do From 94c7409ad5b8812be78a464474f99f41401b32df Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Fri, 10 Feb 2017 14:48:28 -0500 Subject: [PATCH 11/31] drop comments --- src/System/Random/SplitMix/Internal.hs | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index 6c6dac872..420bddf8c 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -136,28 +136,6 @@ instance Monad m => Monad (RandomT m) where (a,_boring) <- ma fseed unRandomT# (f a) nextSeed ---instance Monad m => Monad (RandomT m) where --- (>>=) = \ (RandomT# ma) mf -> --- RandomT# $ \ s -> --- let --- (# splitSeed, nextSeed #) = splitGeneratorSplitMix s --- (# maRes, _boringSeed #) = ma splitSeed --- (# mfRes , resultSeed #) - -{- -there are two models of RandomT m a we could do - -1) s -> (m a , s) - -or - -2) s -> m (a,s) - --- The 'return' function leaves the state unchanged, while @>>=@ uses --- split on the rng state so that the final state of the first computation --- is independent of the second ... -so lets try writing an instance using 1 --} From cfe3c3b21aefbac9538f52b7e3e727ee8a0744af Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Fri, 10 Feb 2017 14:57:43 -0500 Subject: [PATCH 12/31] basic sampler values --- src/System/Random/SplitMix/Internal.hs | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index 420bddf8c..eb290e667 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -8,7 +8,9 @@ module System.Random.SplitMix.Internal( ,nextWord64SplitMix ,SplitMix64(..) ,Random(..) + ,sampleWord64Random ,RandomT(..) + ,sampleWord64RandomT ) where import qualified Data.Bits as DB @@ -109,7 +111,8 @@ instance Monad Random where (# a, _boringSeed #) = ma splitSeed in unRandom# (f a) nextSeed - +sampleWord64Random :: Random Word64 +sampleWord64Random = Random# nextWord64SplitMix newtype RandomT m a = RandomT# { unRandomT# :: (SplitMix64 -> m (a , SplitMix64) ) } @@ -136,6 +139,10 @@ instance Monad m => Monad (RandomT m) where (a,_boring) <- ma fseed unRandomT# (f a) nextSeed +sampleWord64RandomT :: Applicative m => RandomT m Word64 +sampleWord64RandomT = RandomT# $ \ s -> + let (# w, ngen #) = nextWord64SplitMix s + in pure (w, ngen) From 63b9d30faecf077d52acf0bb77dcc5075ec83539 Mon Sep 17 00:00:00 2001 From: Christopher Chalmers Date: Sat, 18 Feb 2017 22:13:25 +0000 Subject: [PATCH 13/31] Add PCG internal module --- random.cabal | 4 +- src/System/Random/PCG/Internal.hs | 160 ++++++++++++++++++++++++++++++ 2 files changed, 163 insertions(+), 1 deletion(-) create mode 100644 src/System/Random/PCG/Internal.hs diff --git a/random.cabal b/random.cabal index 2bc6e1fa6..f186aad29 100644 --- a/random.cabal +++ b/random.cabal @@ -51,7 +51,9 @@ cabal-version: >=1.10 library -- Modules exported by the library. - exposed-modules: System.Random.SplitMix.Internal + exposed-modules: + System.Random.SplitMix.Internal + System.Random.PCG.Internal -- Modules included in this library but not exported. -- other-modules: diff --git a/src/System/Random/PCG/Internal.hs b/src/System/Random/PCG/Internal.hs new file mode 100644 index 000000000..2a5957df4 --- /dev/null +++ b/src/System/Random/PCG/Internal.hs @@ -0,0 +1,160 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE DeriveGeneric #-} +-- | +-- +-- Standard PCG64 Random Number Generator with chosen streams, written in +-- pure haskell. See for details. +-- +module System.Random.PCG.Internal + ( -- * PCG 64 + PCG64 (..) + , seed + , newPCG64 + + -- * Generating random numbers + , next32 + , next64 + , bounded + + -- * Generator utilities + , advancePCG64 + , split + ) where + +import Data.Bits +import Data.Data +import Foreign +import GHC.Generics + +-- | The multiple sequence varient of the pcg random number generator. +data PCG64 = PCG64 + {-# UNPACK #-} !Word64 -- state + {-# UNPACK #-} !Word64 -- inc + deriving (Show, Ord, Eq, Data, Typeable, Generic) + +-- | Create a new generator from two words. +newPCG64 :: Word64 -> Word64 -> PCG64 +newPCG64 = \a b -> + let s = state (PCG64 (a + i) i) + i = (b `shiftL` 1) .|. 1 + in PCG64 s i +{-# INLINE newPCG64 #-} + +-- | Fixed seed. +seed :: PCG64 +seed = PCG64 9600629759793949339 15726070495360670683 +{-# INLINE seed #-} + +-- Generating random numbers ------------------------------------------- + +-- All operations are done via Pair to ensure everything's strict. GHC +-- is good at inlining this so Pair doesn't exist in the generated core. +data Pair = Pair + {-# UNPACK #-} !Word64 -- new state + {-# UNPACK #-} !Word32 -- output + +multiplier :: Word64 +multiplier = 6364136223846793005 +{-# INLINE multiplier #-} + +-- A single step in the generator state +state :: PCG64 -> Word64 +state (PCG64 s inc) = s * multiplier + inc +{-# INLINE state #-} + +-- The random number output +output :: Word64 -> Word32 +output s = + (shifted `unsafeShiftR` rot) .|. (shifted `unsafeShiftL` (negate rot .&. 31)) + where + rot = fromIntegral $ s `shiftR` 59 :: Int + shifted = fromIntegral $ ((s `shiftR` 18) `xor` s) `shiftR` 27 :: Word32 +{-# INLINE output #-} + +-- a new generator state and random number +pair :: PCG64 -> Pair +pair g@(PCG64 s _) = Pair (state g) (output s) +{-# INLINE pair #-} + +-- | Return a random 'Word32' with a generator incremented by one step. +next32 :: PCG64 -> (Word32, PCG64) +next32 = \g@(PCG64 _ inc) -> + let Pair s' r = pair g + in (r, PCG64 s' inc) +{-# INLINE next32 #-} + +-- | Return a random 'Word64' with the generator incremented by two steps. +next64 :: PCG64 -> (Word64, PCG64) +next64 = \(PCG64 s inc) -> + let Pair s' w1 = pair (PCG64 s inc) + Pair s'' w2 = pair (PCG64 s' inc) + in (wordsTo64Bit w1 w2, PCG64 s'' inc) +{-# INLINE next64 #-} + +bounded' :: Word32 -> PCG64 -> Pair +bounded' b (PCG64 s0 inc) = go s0 + where + t = negate b `mod` b + go !s | r >= t = Pair s' (r `mod` b) + | otherwise = go s' + where Pair s' r = pair (PCG64 s inc) +{-# INLINE bounded' #-} + +-- | Generate a uniform 'Word32' less than the bound. If the bound is +-- zero, this throws a divide by zero exception. +bounded :: Word32 -> PCG64 -> (Word32, PCG64) +bounded = \b g@(PCG64 _ inc) -> + let !(Pair s1 w) = bounded' b g + in (w, PCG64 s1 inc) +{-# INLINE bounded #-} + +-- Utilities ----------------------------------------------------------- + +-- | Split the generator @g@ into @(g', g2)@ where @g'@ is @g@ +-- incremented by 4 steps and @g2@ is a new generator with a +-- difference sequence to @g@. +split :: PCG64 -> (PCG64, PCG64) +split (PCG64 s inc) = (PCG64 s4 inc, mk w1 w2 w3 w4) + where + -- This is just something I made up. It passed big crunch and + -- dieharder (by splitting every step) but there's probably a better + -- way. + mk a b c d = newPCG64 (wordsTo64Bit a b) (wordsTo64Bit c d) + Pair s1 w1 = pair (PCG64 s inc) + Pair s2 w2 = pair (PCG64 s1 inc) + Pair s3 w3 = pair (PCG64 s2 inc) + Pair s4 w4 = pair (PCG64 s3 inc) +{-# INLINE split #-} + +advancing + :: Word64 -- amount to advance by + -> Word64 -- state + -> Word64 -- multiplier + -> Word64 -- increment + -> Word64 -- new state +advancing d0 s m0 p0 = go d0 m0 p0 1 0 + where + go !d !cm !cp !am !ap + | d <= 0 = am * s + ap + | odd d = go d' cm' cp' (am * cm) (ap * cm + cp) + | otherwise = go d' cm' cp' am ap + where + cm' = cm * cm + cp' = (cm + 1) * cp + d' = d `div` 2 +{-# INLINE advancing #-} + +-- | Advance a pcg generator @n@ steps. You can give this @-n@ to take +-- the generator back @n@ steps. +advancePCG64 :: Word64 -> PCG64 -> PCG64 +advancePCG64 = \d (PCG64 s inc) -> PCG64 (advancing d s multiplier inc) inc +{-# INLINE advancePCG64 #-} + +-- Misc ---------------------------------------------------------------- + +wordsTo64Bit :: Integral a => Word32 -> Word32 -> a +wordsTo64Bit x y = + fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64) +{-# INLINE wordsTo64Bit #-} + From 6b744c3d31112698718f232b7b3d5986300e9602 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 17 Apr 2017 15:16:33 -0400 Subject: [PATCH 14/31] formmatting tweaks to splitmix code plus enabling wall --- random.cabal | 2 +- src/System/Random/SplitMix/Internal.hs | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/random.cabal b/random.cabal index f186aad29..8e53ba254 100644 --- a/random.cabal +++ b/random.cabal @@ -67,7 +67,7 @@ library -- Directories containing source files. hs-source-dirs: src - ghc-options: -O2 + ghc-options: -O2 -Wall -- 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 index eb290e667..56747211b 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -15,7 +15,7 @@ module System.Random.SplitMix.Internal( import qualified Data.Bits as DB import Data.Bits (xor,(.|.)) -import Data.Word(Word64(..)) +import Data.Word(Word64) import Data.Functor.Identity {-# SPECIALIZE popCount :: Word64 -> Word64 #-} @@ -48,7 +48,6 @@ mix64variant13 = \ w -> xorShiftR 31 $ secondRoundMix64Variant13 $ firstRoundMix firstRoundMix64Variant13 :: Word64 -> Word64 firstRoundMix64Variant13 = \ w -> xorShiftR 30 w * 0xbf58476d1ce4e5b9 - secondRoundMix64Variant13 :: Word64 -> Word64 secondRoundMix64Variant13 = \ w -> xorShiftR 27 w * 0x94d049bb133111eb @@ -76,7 +75,6 @@ data SplitMix64 = SplitMix64 { sm64seed :: {-# UNPACK #-} !Word64 ,sm64Gamma :: {-# UNPACK #-} !Word64 } - advanceSplitMix :: SplitMix64 -> SplitMix64 advanceSplitMix (SplitMix64 sd gamma) = SplitMix64 (sd + gamma) gamma From bbf9330777781ac4a3d49d848343c939ad00c2cb Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sun, 21 May 2017 14:35:08 -0400 Subject: [PATCH 15/31] correct the PCG module name to reflect that its the 32bit word sampler. have a candidate travis file, added entropy to the dependencies for now (spoke with TomMD and we agreed that they should be combined sometime) started work on basic distribution generation tools --- .travis.yml | 103 ++++++++++++++++++++++++++++++ random.cabal | 7 +- src/Data/Distribution/Interval | 4 ++ src/Data/Distribution/Normal.hs | 2 + src/Data/Distribution/Permutation | 0 src/System/Random/PCG/Internal.hs | 82 ++++++++++++------------ 6 files changed, 155 insertions(+), 43 deletions(-) create mode 100644 .travis.yml create mode 100644 src/Data/Distribution/Interval create mode 100644 src/Data/Distribution/Normal.hs create mode 100644 src/Data/Distribution/Permutation diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 000000000..d3e6cac1d --- /dev/null +++ b/.travis.yml @@ -0,0 +1,103 @@ +# This Travis job script has been generated by a script via +# +# make_travis_yml_2.hs 'random.cabal' +# +# For more information, see https://github.com/hvr/multi-ghc-travis +# +language: c +sudo: false + +git: + submodules: false # whether to recursively clone submodules + +cache: + directories: + - $HOME/.cabal/packages + - $HOME/.cabal/store + +before_cache: + - rm -fv $HOME/.cabal/packages/hackage.haskell.org/build-reports.log + # remove files that are regenerated by 'cabal update' + - rm -fv $HOME/.cabal/packages/hackage.haskell.org/00-index.* + - rm -fv $HOME/.cabal/packages/hackage.haskell.org/*.json + - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.cache + - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar + - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar.idx + +matrix: + include: + + + - compiler: "ghc-7.10.3" + env: + # env: TEST=--disable-tests BENCH=--disable-benchmarks + addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-7.10.3], sources: [hvr-ghc]}} + - compiler: "ghc-8.0.1" + env: + # env: TEST=--disable-tests BENCH=--disable-benchmarks + addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.0.1], sources: [hvr-ghc]}} + - compiler: "ghc-8.0.2" + env: + # env: TEST=--disable-tests BENCH=--disable-benchmarks + addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.0.2], sources: [hvr-ghc]}} + - compiler: "ghc-8.2.1" + env: + # env: TEST=--disable-tests BENCH=--disable-benchmarks + addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-8.2.1], sources: [hvr-ghc]}} + - compiler: "ghc-head" + env: + # env: TEST=--disable-tests BENCH=--disable-benchmarks + addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-head], sources: [hvr-ghc]}} + + allow_failures: + - compiler: "ghc-head" + +before_install: + - HC=${CC} + - unset CC + - PATH=/opt/ghc/bin:/opt/ghc-ppa-tools/bin:$PATH + - PKGNAME='random' + +install: + - cabal --version + - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]" + - BENCH=${BENCH---enable-benchmarks} + - TEST=${TEST---enable-tests} + - travis_retry cabal update -v + - sed -i 's/^jobs:/-- jobs:/' ${HOME}/.cabal/config + - rm -fv cabal.project.local + - "echo 'packages: .' > cabal.project" + - rm -f cabal.project.freeze + - cabal new-build -w ${HC} ${TEST} ${BENCH} --dep -j2 + - cabal new-build -w ${HC} --disable-tests --disable-benchmarks --dep -j2 + +# Here starts the actual work to be performed for the package under test; +# any command which exits with a non-zero exit code causes the build to fail. +script: + - if [ -f configure.ac ]; then autoreconf -i; fi + - rm -rf dist/ + - cabal sdist # test that a source-distribution can be generated + - cd dist/ + - SRCTAR=(${PKGNAME}-*.tar.gz) + - SRC_BASENAME="${SRCTAR/%.tar.gz}" + - tar -xvf "./$SRC_BASENAME.tar.gz" + - cd "$SRC_BASENAME/" +## from here on, CWD is inside the extracted source-tarball + - rm -fv cabal.project.local + - "echo 'packages: .' > cabal.project" + # this builds all libraries and executables (without tests/benchmarks) + - rm -f cabal.project.freeze + - cabal new-build -w ${HC} --disable-tests --disable-benchmarks + # this builds all libraries and executables (including tests/benchmarks) + # - rm -rf ./dist-newstyle + - cabal new-build -w ${HC} ${TEST} ${BENCH} + + # there's no 'cabal new-test' yet, so let's emulate for now + - TESTS=( $(awk 'tolower($0) ~ /^test-suite / { print $2 }' *.cabal) ) + - if [ "$TEST" != "--enable-tests" ]; then TESTS=(); fi + - shopt -s globstar; + RC=true; for T in ${TESTS[@]}; do echo "== $T =="; + if dist-newstyle/build/**/$SRC_BASENAME/**/build/$T/$T; then echo "= $T OK ="; + else echo "= $T FAILED ="; RC=false; fi; done; $RC + +# EOF diff --git a/random.cabal b/random.cabal index 8e53ba254..a86215cc6 100644 --- a/random.cabal +++ b/random.cabal @@ -50,7 +50,7 @@ cabal-version: >=1.10 library - -- Modules exported by the library. + -- Modules exported` by the library. exposed-modules: System.Random.SplitMix.Internal System.Random.PCG.Internal @@ -62,7 +62,10 @@ library -- other-extensions: -- Other library packages from which modules are imported. - build-depends: base >=4.8 && <4.10, ghc-prim + build-depends: base >=4.8 && <4.12 + ,ghc-prim + ,entropy == 0.3.* + -- entropy will later be folded into random, probably -- Directories containing source files. hs-source-dirs: src diff --git a/src/Data/Distribution/Interval b/src/Data/Distribution/Interval new file mode 100644 index 000000000..f164b6e9d --- /dev/null +++ b/src/Data/Distribution/Interval @@ -0,0 +1,4 @@ +module Data.Distribution.Interval where + + +unit diff --git a/src/Data/Distribution/Normal.hs b/src/Data/Distribution/Normal.hs new file mode 100644 index 000000000..a7f940696 --- /dev/null +++ b/src/Data/Distribution/Normal.hs @@ -0,0 +1,2 @@ +module Data.Distribution.Normal where + diff --git a/src/Data/Distribution/Permutation b/src/Data/Distribution/Permutation new file mode 100644 index 000000000..e69de29bb diff --git a/src/System/Random/PCG/Internal.hs b/src/System/Random/PCG/Internal.hs index 2a5957df4..e828840f2 100644 --- a/src/System/Random/PCG/Internal.hs +++ b/src/System/Random/PCG/Internal.hs @@ -3,14 +3,14 @@ {-# LANGUAGE DeriveGeneric #-} -- | -- --- Standard PCG64 Random Number Generator with chosen streams, written in +-- Standard PCG32 Random Number Generator with chosen streams, written in -- pure haskell. See for details. -- module System.Random.PCG.Internal - ( -- * PCG 64 - PCG64 (..) + ( -- * PCG 32 + PCG32 (..) , seed - , newPCG64 + , newPCG32 -- * Generating random numbers , next32 @@ -18,7 +18,7 @@ module System.Random.PCG.Internal , bounded -- * Generator utilities - , advancePCG64 + , advancePCG32 , split ) where @@ -28,22 +28,22 @@ import Foreign import GHC.Generics -- | The multiple sequence varient of the pcg random number generator. -data PCG64 = PCG64 +data PCG32 = PCG32 {-# UNPACK #-} !Word64 -- state {-# UNPACK #-} !Word64 -- inc deriving (Show, Ord, Eq, Data, Typeable, Generic) -- | Create a new generator from two words. -newPCG64 :: Word64 -> Word64 -> PCG64 -newPCG64 = \a b -> - let s = state (PCG64 (a + i) i) +newPCG32 :: Word64 -> Word64 -> PCG32 +newPCG32 = \a b -> + let s = state (PCG32 (a + i) i) i = (b `shiftL` 1) .|. 1 - in PCG64 s i -{-# INLINE newPCG64 #-} + in PCG32 s i +{-# INLINE newPCG32 #-} -- | Fixed seed. -seed :: PCG64 -seed = PCG64 9600629759793949339 15726070495360670683 +seed :: PCG32 +seed = PCG32 9600629759793949339 15726070495360670683 {-# INLINE seed #-} -- Generating random numbers ------------------------------------------- @@ -59,8 +59,8 @@ multiplier = 6364136223846793005 {-# INLINE multiplier #-} -- A single step in the generator state -state :: PCG64 -> Word64 -state (PCG64 s inc) = s * multiplier + inc +state :: PCG32 -> Word64 +state (PCG32 s inc) = s * multiplier + inc {-# INLINE state #-} -- The random number output @@ -73,40 +73,40 @@ output s = {-# INLINE output #-} -- a new generator state and random number -pair :: PCG64 -> Pair -pair g@(PCG64 s _) = Pair (state g) (output s) +pair :: PCG32 -> Pair +pair g@(PCG32 s _) = Pair (state g) (output s) {-# INLINE pair #-} -- | Return a random 'Word32' with a generator incremented by one step. -next32 :: PCG64 -> (Word32, PCG64) -next32 = \g@(PCG64 _ inc) -> +next32 :: PCG32 -> (Word32, PCG32) +next32 = \g@(PCG32 _ inc) -> let Pair s' r = pair g - in (r, PCG64 s' inc) + in (r, PCG32 s' inc) {-# INLINE next32 #-} -- | Return a random 'Word64' with the generator incremented by two steps. -next64 :: PCG64 -> (Word64, PCG64) -next64 = \(PCG64 s inc) -> - let Pair s' w1 = pair (PCG64 s inc) - Pair s'' w2 = pair (PCG64 s' inc) - in (wordsTo64Bit w1 w2, PCG64 s'' inc) +next64 :: PCG32 -> (Word64, PCG32) +next64 = \(PCG32 s inc) -> + let Pair s' w1 = pair (PCG32 s inc) + Pair s'' w2 = pair (PCG32 s' inc) + in (wordsTo64Bit w1 w2, PCG32 s'' inc) {-# INLINE next64 #-} -bounded' :: Word32 -> PCG64 -> Pair -bounded' b (PCG64 s0 inc) = go s0 +bounded' :: Word32 -> PCG32 -> Pair +bounded' b (PCG32 s0 inc) = go s0 where t = negate b `mod` b go !s | r >= t = Pair s' (r `mod` b) | otherwise = go s' - where Pair s' r = pair (PCG64 s inc) + where Pair s' r = pair (PCG32 s inc) {-# INLINE bounded' #-} -- | Generate a uniform 'Word32' less than the bound. If the bound is -- zero, this throws a divide by zero exception. -bounded :: Word32 -> PCG64 -> (Word32, PCG64) -bounded = \b g@(PCG64 _ inc) -> +bounded :: Word32 -> PCG32 -> (Word32, PCG32) +bounded = \b g@(PCG32 _ inc) -> let !(Pair s1 w) = bounded' b g - in (w, PCG64 s1 inc) + in (w, PCG32 s1 inc) {-# INLINE bounded #-} -- Utilities ----------------------------------------------------------- @@ -114,17 +114,17 @@ bounded = \b g@(PCG64 _ inc) -> -- | Split the generator @g@ into @(g', g2)@ where @g'@ is @g@ -- incremented by 4 steps and @g2@ is a new generator with a -- difference sequence to @g@. -split :: PCG64 -> (PCG64, PCG64) -split (PCG64 s inc) = (PCG64 s4 inc, mk w1 w2 w3 w4) +split :: PCG32 -> (PCG32, PCG32) +split (PCG32 s inc) = (PCG32 s4 inc, mk w1 w2 w3 w4) where -- This is just something I made up. It passed big crunch and -- dieharder (by splitting every step) but there's probably a better -- way. - mk a b c d = newPCG64 (wordsTo64Bit a b) (wordsTo64Bit c d) - Pair s1 w1 = pair (PCG64 s inc) - Pair s2 w2 = pair (PCG64 s1 inc) - Pair s3 w3 = pair (PCG64 s2 inc) - Pair s4 w4 = pair (PCG64 s3 inc) + mk a b c d = newPCG32 (wordsTo64Bit a b) (wordsTo64Bit c d) + Pair s1 w1 = pair (PCG32 s inc) + Pair s2 w2 = pair (PCG32 s1 inc) + Pair s3 w3 = pair (PCG32 s2 inc) + Pair s4 w4 = pair (PCG32 s3 inc) {-# INLINE split #-} advancing @@ -147,9 +147,9 @@ advancing d0 s m0 p0 = go d0 m0 p0 1 0 -- | Advance a pcg generator @n@ steps. You can give this @-n@ to take -- the generator back @n@ steps. -advancePCG64 :: Word64 -> PCG64 -> PCG64 -advancePCG64 = \d (PCG64 s inc) -> PCG64 (advancing d s multiplier inc) inc -{-# INLINE advancePCG64 #-} +advancePCG32 :: Word64 -> PCG32 -> PCG32 +advancePCG32 = \d (PCG32 s inc) -> PCG32 (advancing d s multiplier inc) inc +{-# INLINE advancePCG32 #-} -- Misc ---------------------------------------------------------------- From f437728f7d75973a779562b1bb2a9b5fc6035bcd Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sun, 21 May 2017 15:02:31 -0400 Subject: [PATCH 16/31] fixed up file names --- src/Data/Distribution/{Interval => Interval.hs} | 0 src/Data/Distribution/Permutation | 0 src/Data/Distribution/Permutation.hs | 2 ++ 3 files changed, 2 insertions(+) rename src/Data/Distribution/{Interval => Interval.hs} (100%) delete mode 100644 src/Data/Distribution/Permutation create mode 100644 src/Data/Distribution/Permutation.hs diff --git a/src/Data/Distribution/Interval b/src/Data/Distribution/Interval.hs similarity index 100% rename from src/Data/Distribution/Interval rename to src/Data/Distribution/Interval.hs diff --git a/src/Data/Distribution/Permutation b/src/Data/Distribution/Permutation deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/Data/Distribution/Permutation.hs b/src/Data/Distribution/Permutation.hs new file mode 100644 index 000000000..dd84b04d6 --- /dev/null +++ b/src/Data/Distribution/Permutation.hs @@ -0,0 +1,2 @@ +module Data.Distribution.Permutation where + From c3b75aa29ac1dd5d376ffab9118d169ff90cdbdb Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sun, 21 May 2017 16:08:11 -0400 Subject: [PATCH 17/31] fix case on CHANGELOG.md --- random.cabal | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/random.cabal b/random.cabal index a86215cc6..06b4dec5e 100644 --- a/random.cabal +++ b/random.cabal @@ -43,7 +43,7 @@ build-type: Simple -- Extra files to be distributed with the package, such as examples or a -- README. -extra-source-files: ChangeLog.md +extra-source-files: CHANGELOG.md -- Constraint on the version of Cabal needed to build this package. cabal-version: >=1.10 From 87dea9042c7336ea3df31aacd10aaa6f21602c36 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Fri, 26 May 2017 11:50:44 -0400 Subject: [PATCH 18/31] provisional fast and really sloppy unit interval generator. For basic uses should be good enough. --- cmmbits/floatsAndBits.cmm | 68 +++++++++++ random.cabal | 19 ++++ src/Data/Distribution/FloatingInterval.hs | 96 ++++++++++++++++ src/Data/Distribution/Interval.hs | 4 - src/Data/Distribution/Normal.hs | 3 + src/Data/Random/Utils.hs | 132 ++++++++++++++++++++++ testCast/WordFloat.hs | 31 +++++ 7 files changed, 349 insertions(+), 4 deletions(-) create mode 100644 cmmbits/floatsAndBits.cmm create mode 100644 src/Data/Distribution/FloatingInterval.hs delete mode 100644 src/Data/Distribution/Interval.hs create mode 100644 src/Data/Random/Utils.hs create mode 100644 testCast/WordFloat.hs diff --git a/cmmbits/floatsAndBits.cmm b/cmmbits/floatsAndBits.cmm new file mode 100644 index 000000000..d610ca3f8 --- /dev/null +++ b/cmmbits/floatsAndBits.cmm @@ -0,0 +1,68 @@ +#include "Cmm.h" +#include "MachDeps.h" + +#if WORD_SIZE_IN_BITS == 64 +#define DOUBLE_SIZE_WDS 1 +#else +#define DOUBLE_SIZE_WDS 2 +#endif + +stg_word64ToDoublezhPrivate(I64 w) +{ + D_ d; + P_ ptr; + + STK_CHK_GEN_N (DOUBLE_SIZE_WDS); + + reserve DOUBLE_SIZE_WDS = ptr { + I64[ptr] = w; + d = D_[ptr]; + } + + return (d); +} + +stg_doubleToWord64zhPrivate(D_ d) +{ + I64 w; + P_ ptr; + + STK_CHK_GEN_N (DOUBLE_SIZE_WDS); + + reserve DOUBLE_SIZE_WDS = ptr { + D_[ptr] = d; + w = I64[ptr]; + } + + return (w); +} + +stg_word32ToFloatzhPrivate(W_ w) +{ + F_ f; + P_ ptr; + + STK_CHK_GEN_N (1); + + reserve 1 = ptr { + I32[ptr] = %lobits32(w); + f = F_[ptr]; + } + + return (f); +} + +stg_floatToWord32zhPrivate(F_ f) +{ + W_ w; + P_ ptr; + + STK_CHK_GEN_N (1); + + reserve 1 = ptr { + F_[ptr] = f; + w = TO_W_(I32[ptr]); + } + + return (w); +} diff --git a/random.cabal b/random.cabal index 06b4dec5e..31944b744 100644 --- a/random.cabal +++ b/random.cabal @@ -54,7 +54,12 @@ library exposed-modules: System.Random.SplitMix.Internal System.Random.PCG.Internal + Data.Distribution.FloatingInterval + Data.Distribution.Normal + Data.Distribution.Permutation + Data.Random.Utils + c-sources: cmmbits/floatsAndBits.cmm -- Modules included in this library but not exported. -- other-modules: @@ -75,3 +80,17 @@ library -- Base language which the package is written in. default-language: Haskell2010 + +test-suite word_and_float + type: exitcode-stdio-1.0 + main-is: WordFloat.hs + hs-source-dirs: + testCast + ghc-options: -Wall + build-depends: + base == 4.* + , random + , tasty == 0.11.* + , tasty-hunit == 0.9.* + other-modules: + default-language: Haskell2010 diff --git a/src/Data/Distribution/FloatingInterval.hs b/src/Data/Distribution/FloatingInterval.hs new file mode 100644 index 000000000..194e60b5c --- /dev/null +++ b/src/Data/Distribution/FloatingInterval.hs @@ -0,0 +1,96 @@ + +{- | This module provides both complete and fast (sloppy) +random samples for the unit interval [+0,1) + + +complete in this context means all representable normalized floats are reachable + +fast means +-} + +{-# LANGUAGE ScopedTypeVariables #-} +module Data.Distribution.FloatingInterval where + +import Data.Bits +import Data.Word +import Data.Random.Utils + +{- +for more info, read the IEEE 2008 Floating point spec or wikipedia. +single precision float is also called binary32, +double precision is also called binary64 + +the greatest negative exponent for double precision normalized floats is -1023, and + 53 bits after the implicit MSB of the significand +for single precision its -126, + and + +in both cases, for normalized form (exponent != 0, theres an implicit leading 1) + + +likewise, the exponent component of a floating point number is "biased" +in the sense of being an unsigned integer, with an additive positive correction factor +to map the actual exponent to its representation. To map from representation to value, +subtract the bias constant + +binary32 bias == 127 +binary64 bias == 1023 +-} + +{- | sampleUnitIntervalDoubleM uniformly samples over the [+0,1) interval of + representable floating point numbers + + TODO: so for now the +-} +sampleUnitIntervalDoubleM :: forall m . Monad m => (m Word64) -> m Double +sampleUnitIntervalDoubleM mword = error "finish this" $ mword {-computeMantissa + where + computeMantissa :: m Double + computeMantissa = do + wd <- mword + leading <- return $ countLeadingZeros wd + if leading == 64 + then computeMoreMantissa 64 + else + computeNormalizedSignificandWith leading wd + + computeNormalizedSignificandWith:: Int -> Word64 -> m Double + computeNormalizedSignificandWith leadingZeros rawWord = + error "finish me" mkUnitDouble leadingZeros rawWord + computeMoreMantissa :: Int -> m Double + computeMoreMantissa = error "finish this too" + --- mkDouble takes the positive version of the negative exponent + --- and the normalized significand (which ellides the leading 1 digit) + mkUnitDouble :: Word64 -> Word64 -> Double + mkUnitDouble negUnBiasedExponent normalSignificand = toIEEE $ undefined (negUnBiasedExponent ) +-} + +{- | sampleUnitIntervalDoubleReallySloppyM, using the same algorithm as in +http://xoroshiro.di.unimi.it/#remarks, which is also used by the rand package +in rust. It has issues, but its super fast. Generates all the representable floats +the correspond to dyadic (binary) rationals of the form k/2^{−53}. Note that +the lowest order bit will 0. Which is why the lowest order bit of the random word +is then xor'd against the corrected unit interval number in this specific + + +extracted docs from the original site: + #include + + (x >> 11) * (1. / (UINT64_C(1) << 53)) +This conversion guarantees that all dyadic rationals of the form k / 2−53 will be equally likely. Note that this conversion prefers the high bits of x, but you can alternatively use the lowest bits. + +An alternative, faster multiplication-free operation is + + #include + + static inline double to_double(uint64_t x) { + const union { uint64_t i; double d; } u = { .i = UINT64_C(0x3FF) << 52 | x >> 12 }; + return u.d - 1.0; + } +The code above cooks up by bit manipulation a real number in the interval [1..2), and then subtracts one to obtain a real number in the interval [0..1). If x is chosen uniformly among 64-bit integers, d is chosen uniformly among dyadic rationals of the form k / 2−52. + -} +sampleUnitIntervalDoubleReallySloppyM :: forall m . Monad m => (m Word64) -> m Double +sampleUnitIntervalDoubleReallySloppyM mword = do + word <- mword + return $ toIEEE $ (\x -> ( x `xor` (1 .&. word) )) $ fromIEEE $ + toIEEE (0x3FF `unsafeShiftL` 52 .|. (word `unsafeShiftR` 12)) - 1 diff --git a/src/Data/Distribution/Interval.hs b/src/Data/Distribution/Interval.hs deleted file mode 100644 index f164b6e9d..000000000 --- a/src/Data/Distribution/Interval.hs +++ /dev/null @@ -1,4 +0,0 @@ -module Data.Distribution.Interval where - - -unit diff --git a/src/Data/Distribution/Normal.hs b/src/Data/Distribution/Normal.hs index a7f940696..de28105a2 100644 --- a/src/Data/Distribution/Normal.hs +++ b/src/Data/Distribution/Normal.hs @@ -1,2 +1,5 @@ module Data.Distribution.Normal where +{- + +-} diff --git a/src/Data/Random/Utils.hs b/src/Data/Random/Utils.hs new file mode 100644 index 000000000..efb64d74d --- /dev/null +++ b/src/Data/Random/Utils.hs @@ -0,0 +1,132 @@ +{-# LANGUAGE Trustworthy #-} +{-# LANGUAGE CPP + , GHCForeignImportPrim + , MagicHash + , UnboxedTuples + , UnliftedFFITypes + #-} +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} + +#include "MachDeps.h" +module Data.Random.Utils( + castWord64ToDouble + ,castDoubleToWord64 + ,castWord32ToFloat + ,castFloatToWord32 + ,CastIEEE(..)) where +import GHC.Word(Word32(..),Word64(..)) +import GHC.Prim (Word#,Float#,Double#) +import GHC.Types +{- +from commit +aa206346e6f12c9f88fdf051185741761ea88fbb +of the ghc git repo, due for inclusion in ghc 8.4 + +this should be move out of random into its own micro package for pre ghc 8.4 compat +-} + + + + +class CastIEEE word ieee | word -> ieee , ieee -> word where + toIEEE :: word -> ieee + fromIEEE :: ieee -> word + +instance CastIEEE Word32 Float where + {-# INLINE toIEEE #-} + {-# INLINE fromIEEE #-} + toIEEE = castWord32ToFloat + fromIEEE = castFloatToWord32 +instance CastIEEE Word64 Double where + {-# INLINE toIEEE #-} + {-# INLINE fromIEEE #-} + toIEEE = castWord64ToDouble + fromIEEE = castDoubleToWord64 + + + + +{- +Note [Casting from integral to floating point types] +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +To implement something like `reinterpret_cast` from C++ to go from a +floating-point type to an integral type one might niavely think that the +following should work: + cast :: Float -> Word32 + cast (F# f#) = W32# (unsafeCoerce# f#) +Unfortunately that is not the case, because all the `unsafeCoerce#` does is tell +the compiler that the types have changed. When one does the above cast and +tries to operate on the resulting `Word32` the code generator will generate code +that performs an integer/word operation on a floating-point register, which +results in a compile error. +The correct way of implementing `reinterpret_cast` to implement a primpop, but +that requires a unique implementation for all supported archetectures. The next +best solution is to write the value from the source register to memory and then +read it from memory into the destination register and the best way to do that +is using CMM. +-} + + + + + + +-- | @'castWord32ToFloat' w@ does a bit-for-bit copy from an integral value +-- to a floating-point value. +-- +-- @since 4.10.0.0 + +{-# INLINE castWord32ToFloat #-} +castWord32ToFloat :: Word32 -> Float +castWord32ToFloat (W32# w#) = F# (stgWord32ToFloat w#) + +foreign import prim "stg_word32ToFloatzhPrivate" + stgWord32ToFloat :: Word# -> Float# + + +-- | @'castFloatToWord32' f@ does a bit-for-bit copy from a floating-point value +-- to an integral value. +-- +-- @since 4.10.0.0 + +{-# INLINE castFloatToWord32 #-} +castFloatToWord32 :: Float -> Word32 +castFloatToWord32 (F# f#) = W32# (stgFloatToWord32 f#) + +foreign import prim "stg_floatToWord32zhPrivate" + stgFloatToWord32 :: Float# -> Word# + + + +-- | @'castWord64ToDouble' w@ does a bit-for-bit copy from an integral value +-- to a floating-point value. +-- +-- @since 4.10.0.0 + +{-# INLINE castWord64ToDouble #-} +castWord64ToDouble :: Word64 -> Double +castWord64ToDouble (W64# w) = D# (stgWord64ToDouble w) + +foreign import prim "stg_word64ToDoublezhPrivate" +#if WORD_SIZE_IN_BITS == 64 + stgWord64ToDouble :: Word# -> Double# +#else + stgWord64ToDouble :: Word64# -> Double# +#endif + + +-- | @'castFloatToWord32' f@ does a bit-for-bit copy from a floating-point value +-- to an integral value. +-- +-- @since 4.10.0.0 + +{-# INLINE castDoubleToWord64 #-} +castDoubleToWord64 :: Double -> Word64 +castDoubleToWord64 (D# d#) = W64# (stgDoubleToWord64 d#) + +foreign import prim "stg_doubleToWord64zhPrivate" +#if WORD_SIZE_IN_BITS == 64 + stgDoubleToWord64 :: Double# -> Word# +#else + stgDoubleToWord64 :: Double# -> Word64# +#endif diff --git a/testCast/WordFloat.hs b/testCast/WordFloat.hs new file mode 100644 index 000000000..b4dea2179 --- /dev/null +++ b/testCast/WordFloat.hs @@ -0,0 +1,31 @@ + +module Main where + + +import Data.Random.Utils +import Data.Word(Word32,Word64) +import Test.Tasty +import Test.Tasty.HUnit + + +main :: IO () +main = defaultMain $ testGroup "zero Rep tests" + [testGroup "single precision" + [testCase "signed zeros representation" $ + assertBool "impossible, binary rep of +/-0 cannot agree" + (castFloatToWord32 (negate 0) /= castFloatToWord32 0) + + ,testCase "zero word32 is +zero" $ + assertBool "word32 zero is +zero, but this failure says nope" + (castWord32ToFloat (0 :: Word32) == (0 :: Float ))] + ,testGroup "double precision zero test " + [testCase "signed zeros" $ + assertBool "zeros agree" + (castDoubleToWord64 (negate 0) /= castDoubleToWord64 0) + + ,testCase "zero word64 is +zero" $ + assertBool "word64 zero is +zero" + (castWord64ToDouble (0 :: Word64) == (0 :: Double )) + ]] + + From a120d4abb55b52332e7a3a2e5115e3301dc1d072 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sat, 27 May 2017 17:07:30 -0400 Subject: [PATCH 19/31] basic polar method rejection sampler for normal distribution --- random.cabal | 1 + src/Data/Distribution/FloatingInterval.hs | 4 +++ src/Data/Distribution/Normal.hs | 33 ++++++++++++++++++++++- 3 files changed, 37 insertions(+), 1 deletion(-) diff --git a/random.cabal b/random.cabal index 31944b744..597865ca0 100644 --- a/random.cabal +++ b/random.cabal @@ -70,6 +70,7 @@ library build-depends: base >=4.8 && <4.12 ,ghc-prim ,entropy == 0.3.* + ,numeric-extras == 0.1.* -- entropy will later be folded into random, probably -- Directories containing source files. diff --git a/src/Data/Distribution/FloatingInterval.hs b/src/Data/Distribution/FloatingInterval.hs index 194e60b5c..b9e8860be 100644 --- a/src/Data/Distribution/FloatingInterval.hs +++ b/src/Data/Distribution/FloatingInterval.hs @@ -74,6 +74,7 @@ is then xor'd against the corrected unit interval number in this specific extracted docs from the original site: +""" #include (x >> 11) * (1. / (UINT64_C(1) << 53)) @@ -88,6 +89,9 @@ An alternative, faster multiplication-free operation is return u.d - 1.0; } The code above cooks up by bit manipulation a real number in the interval [1..2), and then subtracts one to obtain a real number in the interval [0..1). If x is chosen uniformly among 64-bit integers, d is chosen uniformly among dyadic rationals of the form k / 2−52. +""" + + -} sampleUnitIntervalDoubleReallySloppyM :: forall m . Monad m => (m Word64) -> m Double sampleUnitIntervalDoubleReallySloppyM mword = do diff --git a/src/Data/Distribution/Normal.hs b/src/Data/Distribution/Normal.hs index de28105a2..afdb433cc 100644 --- a/src/Data/Distribution/Normal.hs +++ b/src/Data/Distribution/Normal.hs @@ -1,5 +1,36 @@ +{-# LANGUAGE ScopedTypeVariables #-} module Data.Distribution.Normal where -{- +import Numeric.Extras (hypot) +{- | For now this will be using the Marsaglia polar method, +though might be better methods to use in exchange for a teeny bit more complexity. + +NB: tail distribution depends on quality of the unit interval generator + +This implementation only returns one of the two random variates, and thus +only needs to generate real samples from [+0,1) rather than (-1,1) + + +if using the x / 2^-53 style uniform interval, you cant get a result r +such that abs(r) > 15, though thats pretty extreme in the tail distribution -} +unitNormalPolarMethodM :: forall m . Monad m => m Double -> m Bool -> m Double +unitNormalPolarMethodM unitSampler boolSampler = getSample + where + getSample :: m Double + getSample = do + x <- unitSampler + y <- unitSampler + sqrtSumSq <- return $ hypot x y + straightSum <- return $ x*x + y * y + if straightSum > 1 || straightSum == 0 + --- the usual condition is x^2 + y^2 > 1, but the same bound holds for the sqrt thereof + then getSample + else + do + boolS <- boolSampler + signed <- return $ if boolS then 1 else -1 + return $! signed * x * (sqrt ( -2 * log straightSum) / sqrtSumSq) + + From 21ff276a76027d4fb27ad8995e4fa1e5d0dd71d6 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 29 May 2017 17:13:55 -0400 Subject: [PATCH 20/31] rejection sampling for uniform intervals of smallish integers (those representable via a Word64) --- random.cabal | 2 ++ src/Data/Distribution/Integers.hs | 30 ++++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+) create mode 100644 src/Data/Distribution/Integers.hs diff --git a/random.cabal b/random.cabal index 597865ca0..678e561a9 100644 --- a/random.cabal +++ b/random.cabal @@ -57,6 +57,7 @@ library Data.Distribution.FloatingInterval Data.Distribution.Normal Data.Distribution.Permutation + Data.Distribution.Integers Data.Random.Utils c-sources: cmmbits/floatsAndBits.cmm @@ -71,6 +72,7 @@ library ,ghc-prim ,entropy == 0.3.* ,numeric-extras == 0.1.* + ,primitive >= 0.6 -- entropy will later be folded into random, probably -- Directories containing source files. diff --git a/src/Data/Distribution/Integers.hs b/src/Data/Distribution/Integers.hs new file mode 100644 index 000000000..1f5f122e1 --- /dev/null +++ b/src/Data/Distribution/Integers.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE ScopedTypeVariables #-} +module Data.Distribution.Integers where + +import Data.Word(Word64) +import Data.Bits + +--sampleWordRange :: Monad m => (m Word64) -> (Word64,Word64) -> Either + +{- | @'sampleWordRangeSimplified' hi@ samples the closed finite interval @[0,hi]@ -} +sampleWordRangeSimplified :: forall m . Monad m => m Word64 -> Word64 -> m Word64 +sampleWordRangeSimplified mwd upper + | upper + 1 == 0 = mwd + -- full 0 ... 2^64-1 range + | popCount (upper + 1) == 1 = do wd <- mwd ; return (wd `mod` (upper + 1)) + -- power of two range of the form 0 ... 2^k -1, for 0 < k < 64 + | otherwise = rejectionSampler + -- we need to do rejection sampling now! + where + rejectionSampler :: m Word64 + rejectionSampler = do awd <- adjustSampleValue + if awd <= upper then return awd + else rejectionSampler + nextPower2Log :: Int + nextPower2Log = (64 - countLeadingZeros upper ) + adjustSampleValue :: m Word64 + adjustSampleValue = if nextPower2Log == 64 + then mwd + else do wd <- mwd ; return (wd `mod` (bit nextPower2Log)) + + From bb6ec02eb100e29f7dc0504625a9bcbe1846999b Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 29 May 2017 20:55:31 -0400 Subject: [PATCH 21/31] can now uniformly sample any Word64 interval of the form [a,b] --- src/Data/Distribution/Integers.hs | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/Data/Distribution/Integers.hs b/src/Data/Distribution/Integers.hs index 1f5f122e1..a004950f5 100644 --- a/src/Data/Distribution/Integers.hs +++ b/src/Data/Distribution/Integers.hs @@ -4,9 +4,27 @@ module Data.Distribution.Integers where import Data.Word(Word64) import Data.Bits ---sampleWordRange :: Monad m => (m Word64) -> (Word64,Word64) -> Either +{- +this module could probably be generalized to support any @Integral a, FiniteBits a@ +or something like it, though my definition of sampleWordRangeSimplified +assumes mod 2^n wrap around +-} -{- | @'sampleWordRangeSimplified' hi@ samples the closed finite interval @[0,hi]@ -} +{- | @'sampleWordRange' wordSampler (lo,hi)@ will return a uniform sample from the closed interval +@[min lo hi, max lo hi]@ -} +sampleWordRange :: Monad m => m Word64 -> (Word64,Word64) -> m Word64 +sampleWordRange mword (lo,hi) + | lo == hi = return lo + | otherwise = + do wd <- (sampleWordRangeSimplified mword difference) ; return (realLo + wd ) + where + realLo = min lo hi + realHi = max lo hi + difference = realHi - realLo + + + +{- | @'sampleWordRangeSimplified' wordSampler hi@ samples the closed finite interval @[0,hi]@ -} sampleWordRangeSimplified :: forall m . Monad m => m Word64 -> Word64 -> m Word64 sampleWordRangeSimplified mwd upper | upper + 1 == 0 = mwd From 92ad942afd6af23471dfb1f3cebf139a50ea5064 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 29 May 2017 22:39:32 -0400 Subject: [PATCH 22/31] design notes on the permutation algorithm --- src/Data/Distribution/Permutation.hs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/Data/Distribution/Permutation.hs b/src/Data/Distribution/Permutation.hs index dd84b04d6..83a8b48ec 100644 --- a/src/Data/Distribution/Permutation.hs +++ b/src/Data/Distribution/Permutation.hs @@ -1,2 +1,17 @@ module Data.Distribution.Permutation where +import Control.Monad.Primitive +{- | this permutation algorithm is due to knuth. + +to construct a permutation of n symbols, @0... n-1@ + +initialize a size n array A with position @i@ having value @i@ + +(may choose to precompute the sequence of permutations before setting up the array) + +forM [1 .. n-1] \j -> pick a uniform sample s from the interval [j, n-1], +then swap the values at A[j-1] and A[s] + +return the array A +-} +--samplePermutation :: Monad m => From fe65eaa80b9dab60cdfac645afdcaa42a4269cab Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Wed, 31 May 2017 22:58:06 -0400 Subject: [PATCH 23/31] =?UTF-8?q?my=20permutations=20type=20check=20but=20?= =?UTF-8?q?i=20haven=E2=80=99t=20tested=20them=20for=20correctness=20yet?= =?UTF-8?q?=20:)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Data/Distribution/Permutation.hs | 49 ++++++++++++++++++++++++++-- 1 file changed, 47 insertions(+), 2 deletions(-) diff --git a/src/Data/Distribution/Permutation.hs b/src/Data/Distribution/Permutation.hs index 83a8b48ec..3e556aae8 100644 --- a/src/Data/Distribution/Permutation.hs +++ b/src/Data/Distribution/Permutation.hs @@ -1,6 +1,11 @@ +{-# LANGUAGE ScopedTypeVariables, GADTs#-} module Data.Distribution.Permutation where -import Control.Monad.Primitive +import Control.Monad.Primitive as Prim +import Data.Primitive.Array as DPA +import Data.Word(Word32) +import Control.Monad.ST +--import Data.Distribution.Integers {- | this permutation algorithm is due to knuth. to construct a permutation of n symbols, @0... n-1@ @@ -13,5 +18,45 @@ forM [1 .. n-1] \j -> pick a uniform sample s from the interval [j, n-1], then swap the values at A[j-1] and A[s] return the array A + +or to quote the fish-yates shuffle entry on wikipedia + +-- To shuffle an array a of n elements (indices 0..n-1): +for i from 0 to n−2 do + j ← random integer such that i ≤ j < n + exchange a[i] and a[j] + +@`samplePermutation` integerSampler size@ for now is limited to allowing permutations over +no more than 2^32 elements, mostly because if you're wanting larger permutations theres likely better +algorithms available -} ---samplePermutation :: Monad m => +samplePermutation :: forall m . Monad m => ((Word32,Word32)->m Word32) -> Word32 -> m (Array Int) +samplePermutation intervalSample wSize + | wSize == 0 || wSize > 2^(31 :: Int) = error "i'm not letting you do 0 or > 2^31 element permutations" + | otherwise = do + swapList :: [(Int,Int)] <- + mapM (\i -> do jay <- intervalSample (i,wSize - 1) ; return (fromIntegral i,fromIntegral jay ) ) [0 .. wSize - 2 ] + initArray :: Array Int <- return $ fromListN (fromIntegral wSize) [0 .. fromIntegral (wSize - 1)] + return $ runST $ do theMarr <- (thawArray initArray 0 (fromIntegral wSize)) + executeArraySwaps swapList theMarr + +executeArraySwaps :: forall s m . (s~PrimState m,PrimMonad m) => [(Int,Int)] + -> MutableArray s Int -> m (Array Int) +executeArraySwaps [] _marr = error "you really shouldn't be invoking executeArraySwaps on empty things" +executeArraySwaps ls@((a,_):_) marr + | a /= 0 = error "the swap sequence list for executeArraySwaps doesn't start with a swap with zero, this is a bug!" + | otherwise = do swapSpots 0 ls ; unsafeFreezeArray marr + where + arrayLength :: Int + arrayLength = sizeofMutableArray marr + swapSpots :: Int -> [(Int,Int)] -> m () + swapSpots ix [] | ix >= (arrayLength - 2) = return () + | otherwise = error "the swap list for executeArraySwaps is shorter than the array length" + swapSpots ix _ | ix >= arrayLength = error "can't swap permutations beyond the array size in executeArraySwaps" + swapSpots ix ((from,to):rest) + | ix /= from = error "bad coordinate mismatch " + | otherwise = do aval <- readArray marr from ; bval <- readArray marr to + writeArray marr from bval ; writeArray marr to aval + swapSpots (ix +1) rest + + From 615434ec425322ff2cab7d33e1c6211a91691972 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Wed, 31 May 2017 22:58:26 -0400 Subject: [PATCH 24/31] lots of miscellaneous cleanups --- random.cabal | 9 ++++++++- readme.md | 5 +++++ src/Data/Distribution/FloatingInterval.hs | 3 ++- src/Data/Distribution/Integers.hs | 4 +++- src/Data/Distribution/Normal.hs | 4 ++-- 5 files changed, 20 insertions(+), 5 deletions(-) create mode 100644 readme.md diff --git a/random.cabal b/random.cabal index 678e561a9..027cbf40f 100644 --- a/random.cabal +++ b/random.cabal @@ -43,7 +43,12 @@ build-type: Simple -- Extra files to be distributed with the package, such as examples or a -- README. -extra-source-files: CHANGELOG.md +extra-source-files: + CHANGELOG.md + .travis.yml + readme.md + + -- Constraint on the version of Cabal needed to build this package. cabal-version: >=1.10 @@ -73,6 +78,8 @@ library ,entropy == 0.3.* ,numeric-extras == 0.1.* ,primitive >= 0.6 + ,transformers >= 0.2 + ,transformers-compat >= 0.3 -- entropy will later be folded into random, probably -- Directories containing source files. diff --git a/readme.md b/readme.md new file mode 100644 index 000000000..e87917c42 --- /dev/null +++ b/readme.md @@ -0,0 +1,5 @@ +https://travis-ci.org/cartazio/random.svg?branch=master + +# Random + +a random number generation and sampling library for haskell diff --git a/src/Data/Distribution/FloatingInterval.hs b/src/Data/Distribution/FloatingInterval.hs index b9e8860be..62cda7514 100644 --- a/src/Data/Distribution/FloatingInterval.hs +++ b/src/Data/Distribution/FloatingInterval.hs @@ -37,6 +37,7 @@ binary32 bias == 127 binary64 bias == 1023 -} + {- | sampleUnitIntervalDoubleM uniformly samples over the [+0,1) interval of representable floating point numbers @@ -78,7 +79,7 @@ extracted docs from the original site: #include (x >> 11) * (1. / (UINT64_C(1) << 53)) -This conversion guarantees that all dyadic rationals of the form k / 2−53 will be equally likely. Note that this conversion prefers the high bits of x, but you can alternatively use the lowest bits. +This conversion guarantees that all dyadic rationals of the form k / 2^53 will be equally likely. Note that this conversion prefers the high bits of x, but you can alternatively use the lowest bits. An alternative, faster multiplication-free operation is diff --git a/src/Data/Distribution/Integers.hs b/src/Data/Distribution/Integers.hs index a004950f5..e0d4959e5 100644 --- a/src/Data/Distribution/Integers.hs +++ b/src/Data/Distribution/Integers.hs @@ -11,7 +11,9 @@ assumes mod 2^n wrap around -} {- | @'sampleWordRange' wordSampler (lo,hi)@ will return a uniform sample from the closed interval -@[min lo hi, max lo hi]@ -} +@[min lo hi, max lo hi]@ +maybe should throw error instead, not sure :) +-} sampleWordRange :: Monad m => m Word64 -> (Word64,Word64) -> m Word64 sampleWordRange mword (lo,hi) | lo == hi = return lo diff --git a/src/Data/Distribution/Normal.hs b/src/Data/Distribution/Normal.hs index afdb433cc..4c25e38f5 100644 --- a/src/Data/Distribution/Normal.hs +++ b/src/Data/Distribution/Normal.hs @@ -22,8 +22,8 @@ unitNormalPolarMethodM unitSampler boolSampler = getSample getSample = do x <- unitSampler y <- unitSampler - sqrtSumSq <- return $ hypot x y - straightSum <- return $ x*x + y * y + sqrtSumSq <- return $ hypot x y + straightSum <- return $ x*x + y*y if straightSum > 1 || straightSum == 0 --- the usual condition is x^2 + y^2 > 1, but the same bound holds for the sqrt thereof then getSample From 253eab23081f1603ad9c4aa53209cb59d4a01778 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Sat, 3 Jun 2017 16:07:22 -0400 Subject: [PATCH 25/31] add vector dependency for now and clean up the permutation code and its representation shifting to unboxed vector for the permutations to improve locality and memory usage --- random.cabal | 2 ++ src/Data/Distribution/Permutation.hs | 44 +++++++++++++++++----------- 2 files changed, 29 insertions(+), 17 deletions(-) diff --git a/random.cabal b/random.cabal index 027cbf40f..9b81aa653 100644 --- a/random.cabal +++ b/random.cabal @@ -80,6 +80,8 @@ library ,primitive >= 0.6 ,transformers >= 0.2 ,transformers-compat >= 0.3 + ,vector >= 0.11 && < 0.14 + -- vector 0.13 wont likely break stuff i care about -- entropy will later be folded into random, probably -- Directories containing source files. diff --git a/src/Data/Distribution/Permutation.hs b/src/Data/Distribution/Permutation.hs index 3e556aae8..929b78011 100644 --- a/src/Data/Distribution/Permutation.hs +++ b/src/Data/Distribution/Permutation.hs @@ -2,9 +2,14 @@ module Data.Distribution.Permutation where import Control.Monad.Primitive as Prim -import Data.Primitive.Array as DPA +--import Data.Primitive.Array as DPA import Data.Word(Word32) -import Control.Monad.ST +import Data.Int(Int32) +import Control.Monad.ST(runST) +import Control.Monad(forM,forM_) +import Data.Vector.Unboxed.Mutable as DVM +import qualified Data.Vector.Unboxed as DV + --import Data.Distribution.Integers {- | this permutation algorithm is due to knuth. @@ -30,33 +35,38 @@ for i from 0 to n−2 do no more than 2^32 elements, mostly because if you're wanting larger permutations theres likely better algorithms available -} -samplePermutation :: forall m . Monad m => ((Word32,Word32)->m Word32) -> Word32 -> m (Array Int) +samplePermutation :: forall m . (Monad m )=> ((Word32,Word32)->m Word32) -> Word32 -> m (DV.Vector Int32) samplePermutation intervalSample wSize | wSize == 0 || wSize > 2^(31 :: Int) = error "i'm not letting you do 0 or > 2^31 element permutations" | otherwise = do - swapList :: [(Int,Int)] <- - mapM (\i -> do jay <- intervalSample (i,wSize - 1) ; return (fromIntegral i,fromIntegral jay ) ) [0 .. wSize - 2 ] - initArray :: Array Int <- return $ fromListN (fromIntegral wSize) [0 .. fromIntegral (wSize - 1)] - return $ runST $ do theMarr <- (thawArray initArray 0 (fromIntegral wSize)) - executeArraySwaps swapList theMarr + swapList :: [(Int,Int)] <- forM [0 .. wSize - 2 ] + (\i -> do jay <- intervalSample (i,wSize - 1) ; + return (fromIntegral i,fromIntegral jay ) ) + + return $ runST $ + do vecM <- DVM.unsafeNew (fromIntegral wSize) + forM_ [0 :: Int .. fromIntegral wSize - 1 ] + (\ i -> DVM.write vecM i (fromIntegral i :: Int32)) + executeArraySwaps swapList vecM executeArraySwaps :: forall s m . (s~PrimState m,PrimMonad m) => [(Int,Int)] - -> MutableArray s Int -> m (Array Int) + -> DVM.MVector s Int32 -> m (DV.Vector Int32) executeArraySwaps [] _marr = error "you really shouldn't be invoking executeArraySwaps on empty things" executeArraySwaps ls@((a,_):_) marr - | a /= 0 = error "the swap sequence list for executeArraySwaps doesn't start with a swap with zero, this is a bug!" - | otherwise = do swapSpots 0 ls ; unsafeFreezeArray marr + | a /= 0 = error "the swap sequence list for executeArraySwaps doesn't start with a swap with zero" + | otherwise = do swapSpots 0 ls ; DV.unsafeFreeze marr where arrayLength :: Int - arrayLength = sizeofMutableArray marr + arrayLength = DVM.length marr swapSpots :: Int -> [(Int,Int)] -> m () - swapSpots ix [] | ix >= (arrayLength - 2) = return () - | otherwise = error "the swap list for executeArraySwaps is shorter than the array length" - swapSpots ix _ | ix >= arrayLength = error "can't swap permutations beyond the array size in executeArraySwaps" + swapSpots ix [] + | ix >= (arrayLength - 2) = return () + | otherwise = error "the swap list for executeArraySwaps is shorter than the array length" + swapSpots ix _ + | ix >= (arrayLength - 1 ) = error "can't swap permutations beyond the array size in executeArraySwaps" swapSpots ix ((from,to):rest) | ix /= from = error "bad coordinate mismatch " - | otherwise = do aval <- readArray marr from ; bval <- readArray marr to - writeArray marr from bval ; writeArray marr to aval + | otherwise = do DVM.swap marr from to swapSpots (ix +1) rest From c9aea5b01f6997c5fe0e39834bfa6b54f22d2c15 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 6 Nov 2017 20:32:00 -0500 Subject: [PATCH 26/31] some formatting cleanup and current task list --- .gitignore | 1 + TODO | 15 +++++++++++++++ src/Data/Distribution/Normal.hs | 2 +- src/Data/Distribution/Permutation.hs | 4 ++-- src/System/Random/SplitMix/Internal.hs | 6 ++++-- 5 files changed, 23 insertions(+), 5 deletions(-) create mode 100644 .gitignore create mode 100644 TODO diff --git a/.gitignore b/.gitignore new file mode 100644 index 000000000..c70beef85 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +cabal.project.local diff --git a/TODO b/TODO new file mode 100644 index 000000000..ef0dea00b --- /dev/null +++ b/TODO @@ -0,0 +1,15 @@ +Random TODO + +* have fast gen float / double siblings +* modern travis CI setup +* moving the c code back into haskell +* criterion benchmarks to compare +* test suite to compare impls +* test suite to validate determination of stream given several seeds +* migrating the bad seed mitigation stuff from the Java sibling +* emailing GUY!!!! +* doc updates +* considering UX changes +* RandomT?? + * Should monadic bind do a split (a la quicheck ) +* cleanup big crush stuff and have that as an ancillary library \ No newline at end of file diff --git a/src/Data/Distribution/Normal.hs b/src/Data/Distribution/Normal.hs index 4c25e38f5..910db8cdc 100644 --- a/src/Data/Distribution/Normal.hs +++ b/src/Data/Distribution/Normal.hs @@ -24,7 +24,7 @@ unitNormalPolarMethodM unitSampler boolSampler = getSample y <- unitSampler sqrtSumSq <- return $ hypot x y straightSum <- return $ x*x + y*y - if straightSum > 1 || straightSum == 0 + if straightSum >= 1 || straightSum == 0 --- the usual condition is x^2 + y^2 > 1, but the same bound holds for the sqrt thereof then getSample else diff --git a/src/Data/Distribution/Permutation.hs b/src/Data/Distribution/Permutation.hs index 929b78011..9bd642638 100644 --- a/src/Data/Distribution/Permutation.hs +++ b/src/Data/Distribution/Permutation.hs @@ -24,7 +24,7 @@ then swap the values at A[j-1] and A[s] return the array A -or to quote the fish-yates shuffle entry on wikipedia +or to quote the fisher-yates shuffle entry on wikipedia -- To shuffle an array a of n elements (indices 0..n-1): for i from 0 to n−2 do @@ -35,7 +35,7 @@ for i from 0 to n−2 do no more than 2^32 elements, mostly because if you're wanting larger permutations theres likely better algorithms available -} -samplePermutation :: forall m . (Monad m )=> ((Word32,Word32)->m Word32) -> Word32 -> m (DV.Vector Int32) +samplePermutation :: forall m . (Monad m) => ((Word32,Word32)->m Word32) -> Word32 -> m (DV.Vector Int32) samplePermutation intervalSample wSize | wSize == 0 || wSize > 2^(31 :: Int) = error "i'm not letting you do 0 or > 2^31 element permutations" | otherwise = do diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index 56747211b..14eb41b82 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -139,10 +139,12 @@ instance Monad m => Monad (RandomT m) where sampleWord64RandomT :: Applicative m => RandomT m Word64 sampleWord64RandomT = RandomT# $ \ s -> - let (# w, ngen #) = nextWord64SplitMix s + let (# !w, !ngen #) = nextWord64SplitMix s in pure (w, ngen) - +--instance PrimMonad m => PrimMonad (RandomT m) where +-- primitive = \ m -> +-- {-# INLINE #-} nextWord64SplitMix :: SplitMix64 -> (# Word64 , SplitMix64 #) nextWord64SplitMix gen = mixedRes `seq` (# mixedRes , newgen #) From c85b0d007d81aa9ec122155be361dc26e3f45e58 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 25 Dec 2017 14:19:15 -0500 Subject: [PATCH 27/31] making the mangling safe for ghc 8.4 and above --- cmmbits/floatsAndBits.cmm | 8 ++++---- src/Data/Random/Utils.hs | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/cmmbits/floatsAndBits.cmm b/cmmbits/floatsAndBits.cmm index d610ca3f8..acfa456be 100644 --- a/cmmbits/floatsAndBits.cmm +++ b/cmmbits/floatsAndBits.cmm @@ -7,7 +7,7 @@ #define DOUBLE_SIZE_WDS 2 #endif -stg_word64ToDoublezhPrivate(I64 w) +cts_random_stg_word64ToDoublezhPrivate(I64 w) { D_ d; P_ ptr; @@ -22,7 +22,7 @@ stg_word64ToDoublezhPrivate(I64 w) return (d); } -stg_doubleToWord64zhPrivate(D_ d) +cts_random_stg_doubleToWord64zhPrivate(D_ d) { I64 w; P_ ptr; @@ -37,7 +37,7 @@ stg_doubleToWord64zhPrivate(D_ d) return (w); } -stg_word32ToFloatzhPrivate(W_ w) +cts_random_stg_word32ToFloatzhPrivate(W_ w) { F_ f; P_ ptr; @@ -52,7 +52,7 @@ stg_word32ToFloatzhPrivate(W_ w) return (f); } -stg_floatToWord32zhPrivate(F_ f) +cts_random_stg_floatToWord32zhPrivate(F_ f) { W_ w; P_ ptr; diff --git a/src/Data/Random/Utils.hs b/src/Data/Random/Utils.hs index efb64d74d..560b333a5 100644 --- a/src/Data/Random/Utils.hs +++ b/src/Data/Random/Utils.hs @@ -80,7 +80,7 @@ is using CMM. castWord32ToFloat :: Word32 -> Float castWord32ToFloat (W32# w#) = F# (stgWord32ToFloat w#) -foreign import prim "stg_word32ToFloatzhPrivate" +foreign import prim "cts_random_stg_word32ToFloatzhPrivate" stgWord32ToFloat :: Word# -> Float# @@ -93,7 +93,7 @@ foreign import prim "stg_word32ToFloatzhPrivate" castFloatToWord32 :: Float -> Word32 castFloatToWord32 (F# f#) = W32# (stgFloatToWord32 f#) -foreign import prim "stg_floatToWord32zhPrivate" +foreign import prim "cts_random_stg_floatToWord32zhPrivate" stgFloatToWord32 :: Float# -> Word# @@ -107,7 +107,7 @@ foreign import prim "stg_floatToWord32zhPrivate" castWord64ToDouble :: Word64 -> Double castWord64ToDouble (W64# w) = D# (stgWord64ToDouble w) -foreign import prim "stg_word64ToDoublezhPrivate" +foreign import prim "cts_random_stg_word64ToDoublezhPrivate" #if WORD_SIZE_IN_BITS == 64 stgWord64ToDouble :: Word# -> Double# #else @@ -124,7 +124,7 @@ foreign import prim "stg_word64ToDoublezhPrivate" castDoubleToWord64 :: Double -> Word64 castDoubleToWord64 (D# d#) = W64# (stgDoubleToWord64 d#) -foreign import prim "stg_doubleToWord64zhPrivate" +foreign import prim "cts_random_stg_doubleToWord64zhPrivate" #if WORD_SIZE_IN_BITS == 64 stgDoubleToWord64 :: Double# -> Word# #else From de51f60924d0a04b3c21bd7b778673925aad36bf Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 25 Dec 2017 20:12:23 -0500 Subject: [PATCH 28/31] setting up the skeleton for a 1.3 mostly compatible release --- random.cabal | 6 +- src/Data/Random/Utils.hs | 12 +- src/System/Random.hs | 525 ++++++++++++++++++ src/System/Random/{PCG => PCG32}/Internal.hs | 2 +- src/System/Random/SplitMix/Internal.hs | 41 +- .../Random/SplitMix/Internal/Splitting.hs | 2 + 6 files changed, 568 insertions(+), 20 deletions(-) create mode 100644 src/System/Random.hs rename src/System/Random/{PCG => PCG32}/Internal.hs (99%) create mode 100644 src/System/Random/SplitMix/Internal/Splitting.hs diff --git a/random.cabal b/random.cabal index 9b81aa653..2d4f18d21 100644 --- a/random.cabal +++ b/random.cabal @@ -10,7 +10,7 @@ name: random -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 2.0.0.0 +version: 1.3.0.0 -- A short (one-line) description of the package. synopsis: Random number generation library for haskell @@ -57,8 +57,10 @@ cabal-version: >=1.10 library -- Modules exported` by the library. exposed-modules: + System.Random System.Random.SplitMix.Internal - System.Random.PCG.Internal + System.Random.SplitMix.Internal.Splitting + System.Random.PCG32.Internal Data.Distribution.FloatingInterval Data.Distribution.Normal Data.Distribution.Permutation diff --git a/src/Data/Random/Utils.hs b/src/Data/Random/Utils.hs index 560b333a5..4bd2ec1ad 100644 --- a/src/Data/Random/Utils.hs +++ b/src/Data/Random/Utils.hs @@ -13,7 +13,7 @@ module Data.Random.Utils( ,castDoubleToWord64 ,castWord32ToFloat ,castFloatToWord32 - ,CastIEEE(..)) where + ,RandomCastIEEE(..)) where import GHC.Word(Word32(..),Word64(..)) import GHC.Prim (Word#,Float#,Double#) import GHC.Types @@ -23,21 +23,25 @@ aa206346e6f12c9f88fdf051185741761ea88fbb of the ghc git repo, due for inclusion in ghc 8.4 this should be move out of random into its own micro package for pre ghc 8.4 compat +with conversion facilities in ghc >= 8.4 + +this copy has name mangling at the CMM layer for happy linking +plus Random prefixing the class so it should be low headache -} -class CastIEEE word ieee | word -> ieee , ieee -> word where +class RandomCastIEEE word ieee | word -> ieee , ieee -> word where toIEEE :: word -> ieee fromIEEE :: ieee -> word -instance CastIEEE Word32 Float where +instance RandomCastIEEE Word32 Float where {-# INLINE toIEEE #-} {-# INLINE fromIEEE #-} toIEEE = castWord32ToFloat fromIEEE = castFloatToWord32 -instance CastIEEE Word64 Double where +instance RandomCastIEEE Word64 Double where {-# INLINE toIEEE #-} {-# INLINE fromIEEE #-} toIEEE = castWord64ToDouble diff --git a/src/System/Random.hs b/src/System/Random.hs new file mode 100644 index 000000000..f6bde7f35 --- /dev/null +++ b/src/System/Random.hs @@ -0,0 +1,525 @@ +{-# LANGUAGE CPP #-} +#if __GLASGOW_HASKELL__ >= 701 +{-# LANGUAGE Trustworthy #-} +#endif + +----------------------------------------------------------------------------- +-- | +-- Module : System.Random +-- Copyright : (c) The University of Glasgow 2001 +-- License : BSD-style (see the file LICENSE in the 'random' repository) +-- +-- Maintainer : libraries@haskell.org +-- Stability : stable +-- Portability : portable +-- +-- This library deals with the common task of pseudo-random number +-- generation. The library makes it possible to generate repeatable +-- results, by starting with a specified initial random number generator, +-- or to get different results on each run by using the system-initialised +-- generator or by supplying a seed from some other source. +-- +-- The library is split into two layers: +-- +-- * A core /random number generator/ provides a supply of bits. +-- The class 'RandomGen' provides a common interface to such generators. +-- The library provides one instance of 'RandomGen', the abstract +-- data type 'StdGen'. Programmers may, of course, supply their own +-- instances of 'RandomGen'. +-- +-- * The class 'Random' provides a way to extract values of a particular +-- type from a random number generator. For example, the 'Float' +-- 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. +-- +----------------------------------------------------------------------------- + +#include "MachDeps.h" + +module System.Random + ( + + -- $intro + + -- * Random number generators + + + RandomGen(next, genRange) + , SplittableGen(split) + + -- ** Standard random number generators + --, StdGen + --, mkStdGen + + -- ** The global random number generator + + -- $globalrng + + --, getStdRandom + --, getStdGen + --, setStdGen + --, newStdGen + + -- * Random values of various types + , Random ( random, randomR, + randoms, randomRs, + randomIO, randomRIO ) + + -- * References + -- $references + + ) where + +import Prelude + +import Data.Bits +import Data.Int +import Data.Word +import Foreign.C.Types + + + +--import Data.Ratio ( numerator, denominator ) + +--import Data.Char ( isSpace, chr, ord ) +--import System.IO.Unsafe ( unsafePerformIO ) +--import Data.IORef ( IORef, newIORef, readIORef, writeIORef ) +#if MIN_VERSION_base (4,6,0) +--import Data.IORef ( atomicModifyIORef' ) +#else +import Data.IORef ( atomicModifyIORef ) +#endif +--import Numeric ( readDec ) + +#ifdef __GLASGOW_HASKELL__ +import GHC.Exts ( build ) +#else +-- | A dummy variant of build without fusion. +{-# INLINE build #-} +build :: ((a -> [a] -> [a]) -> [a] -> [a]) -> [a] +build g = g (:) [] +#endif + +#if !MIN_VERSION_base (4,6,0) +atomicModifyIORef' :: IORef a -> (a -> (a,b)) -> IO b +atomicModifyIORef' ref f = do + b <- atomicModifyIORef ref + (\x -> let (a, b) = f x + in (a, a `seq` b)) + b `seq` return b +#endif + + + +-- | The class 'RandomGen' provides a common interface to random number +-- generators. +-- +-- Minimal complete definition: 'next'. + +class RandomGen g where + + -- |The 'next' operation returns an 'Int' that is uniformly distributed + -- in the range returned by 'genRange' (including both end points), + -- and a new generator. + next :: g -> (Word64, g) + + -- |The 'genRange' operation yields the range of values returned by + -- the generator. + -- + -- It is required that: + -- + -- * If @(a,b) = 'genRange' g@, then @a < b@. + -- + -- * 'genRange' always returns a pair of defined 'Int's. + -- + -- The second condition ensures that 'genRange' cannot examine its + -- argument, and hence the value it returns can be determined only by the + -- instance of 'RandomGen'. That in turn allows an implementation to make + -- a single call to 'genRange' to establish a generator's range, without + -- being concerned that the generator returned by (say) 'next' might have + -- a different range to the generator passed to 'next'. + -- + -- The default definition spans the full range of 'Int'. + genRange :: g -> (Word64,Word64) + + -- default method + genRange _ = (minBound, maxBound) + + +-- | The class 'SplittableGen' proivides a way to specify a random number +-- generator that can be split into two new generators. +class SplittableGen g where + + -- |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). + 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@. + +In addition, 'reads' may be used to map an arbitrary string (not necessarily one +produced by 'show') onto a value of type 'StdGen'. In general, the 'Read' +instance of 'StdGen' has the following properties: + +* It guarantees to succeed on any string. + +* 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 +{- +instance RandomGen StdGen where + next = stdNext + genRange _ = stdRange + + +instance SplittableGen StdGen where + 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. +-} +stdFromString :: String -> (StdGen, String) +stdFromString s = (mkStdGen num, rest) + where (cs, rest) = splitAt 6 s + num = foldl (\a x -> x + 3 * a) 1 (map ord cs) + + + | +The function 'mkStdGen' provides an alternative way of producing an initial +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 + +{- +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 + +createStdGen :: Integer -> StdGen +createStdGen s = mkStdGen32 $ fromIntegral s +-} + +{- | +With a source of random number supply in hand, the 'Random' class allows the +programmer to extract random values of a variety of types. + +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. + randomR :: RandomGen g => (a,a) -> g -> (a,g) + + -- | The same as 'randomR', but using a default range determined by the type: + -- + -- * For bounded types (instances of 'Bounded', such as 'Char'), + -- the range is normally the whole type. + -- + -- * For fractional types, the range is normally the semi-closed interval + -- @[0,1)@. + -- + -- * For 'Integer', the range is (arbitrarily) the range of 'Int'. + random :: RandomGen g => g -> (a, g) + + -- | Plural variant of 'randomR', producing an infinite list of + -- random values instead of returning a new generator. + {-# INLINE randomRs #-} + randomRs :: RandomGen g => (a,a) -> g -> [a] + randomRs ival g = build (\cons _nil -> buildRandoms cons (randomR ival) g) + + -- | Plural variant of 'random', producing an infinite list of + -- random values instead of returning a new generator. + {-# INLINE randoms #-} + randoms :: RandomGen g => g -> [a] + randoms g = build (\cons _nil -> buildRandoms cons random g) + + -- | A variant of 'randomR' that uses the global random number generator + -- (see "System.Random#globalrng"). + randomRIO :: (a,a) -> IO a + --randomRIO range = getStdRandom (randomR range) + + -- | A variant of 'random' that uses the global random number generator + -- (see "System.Random#globalrng"). + randomIO :: IO a + --randomIO = getStdRandom random + +-- | Produce an infinite list-equivalent of random values. +{-# INLINE buildRandoms #-} +buildRandoms :: RandomGen g + => (a -> as -> as) -- ^ E.g. '(:)' but subject to fusion + -> (g -> (a,g)) -- ^ E.g. 'random' + -> g -- ^ A 'RandomGen' instance + -> as +buildRandoms cons rand = go + where + -- The seq fixes part of #4218 and also makes fused Core simpler. + go g = x `seq` (x `cons` go g') where (x,g') = rand g + + +instance Random Integer where + randomR ival g = randomIvalInteger ival g + random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g + +instance Random Int where randomR = randomIvalIntegral; random = randomBounded +instance Random Int8 where randomR = randomIvalIntegral; random = randomBounded +instance Random Int16 where randomR = randomIvalIntegral; random = randomBounded +instance Random Int32 where randomR = randomIvalIntegral; random = randomBounded +instance Random Int64 where randomR = randomIvalIntegral; random = randomBounded + + + +instance Random Word where randomR = randomIvalIntegral; random = randomBounded +instance Random Word8 where randomR = randomIvalIntegral; random = randomBounded +instance Random Word16 where randomR = randomIvalIntegral; random = randomBounded +instance Random Word32 where randomR = randomIvalIntegral; random = randomBounded +instance Random Word64 where randomR = randomIvalIntegral; random = randomBounded + +instance Random CChar where randomR = randomIvalIntegral; random = randomBounded +instance Random CSChar where randomR = randomIvalIntegral; random = randomBounded +instance Random CUChar where randomR = randomIvalIntegral; random = randomBounded +instance Random CShort where randomR = randomIvalIntegral; random = randomBounded +instance Random CUShort where randomR = randomIvalIntegral; random = randomBounded +instance Random CInt where randomR = randomIvalIntegral; random = randomBounded +instance Random CUInt where randomR = randomIvalIntegral; random = randomBounded +instance Random CLong where randomR = randomIvalIntegral; random = randomBounded +instance Random CULong where randomR = randomIvalIntegral; random = randomBounded +instance Random CPtrdiff where randomR = randomIvalIntegral; random = randomBounded +instance Random CSize where randomR = randomIvalIntegral; random = randomBounded +instance Random CWchar where randomR = randomIvalIntegral; random = randomBounded +instance Random CSigAtomic where randomR = randomIvalIntegral; random = randomBounded +instance Random CLLong where randomR = randomIvalIntegral; random = randomBounded +instance Random CULLong where randomR = randomIvalIntegral; random = randomBounded +instance Random CIntPtr where randomR = randomIvalIntegral; random = randomBounded +instance Random CUIntPtr where randomR = randomIvalIntegral; random = randomBounded +instance Random CIntMax where randomR = randomIvalIntegral; random = randomBounded +instance Random CUIntMax where randomR = randomIvalIntegral; random = randomBounded + +instance Random Char where + --randomR (a,b) g = + -- case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of + -- (x,g') -> (chr x, g') + --random g = randomR (minBound,maxBound) g + +instance Random Bool where + randomR (a,b) g = + case (randomIvalInteger (bool2Int a, bool2Int b) g) of + (x, g') -> (int2Bool x, g') + where + bool2Int :: Bool -> Integer + bool2Int False = 0 + bool2Int True = 1 + + int2Bool :: Int -> Bool + int2Bool 0 = False + int2Bool _ = True + + random g = randomR (minBound,maxBound) g + +{-# INLINE randomRFloating #-} +randomRFloating :: (Fractional a, Num a, Ord a, Random a, RandomGen g) => (a, a) -> g -> (a, g) +randomRFloating (l,h) g + | l>h = randomRFloating (h,l) g + | otherwise = let (coef,g') = random g in + (2.0 * (0.5*l + coef * (0.5*h - 0.5*l)), g') -- avoid overflow + +instance Random Double where + randomR = randomRFloating + random rng = + case random rng of + (x,rng') -> + -- We use 53 bits of randomness corresponding to the 53 bit significand: + ((fromIntegral (mask53 .&. (x::Int64)) :: Double) + / fromIntegral twoto53, rng') + where + twoto53 = (2::Int64) ^ (53::Int64) + mask53 = twoto53 - 1 + +instance Random Float where + randomR = randomRFloating + random rng = + -- TODO: Faster to just use 'next' IF it generates enough bits of randomness. + case random rng of + (x,rng') -> + -- We use 24 bits of randomness corresponding to the 24 bit significand: + ((fromIntegral (mask24 .&. (x::Int32)) :: Float) + / fromIntegral twoto24, rng') + -- Note, encodeFloat is another option, but I'm not seeing slightly + -- worse performance with the following [2011.06.25]: +-- (encodeFloat rand (-24), rng') + where + mask24 = twoto24 - 1 + twoto24 = (2::Int32) ^ (24::Int32) + +-- CFloat/CDouble are basically the same as a Float/Double: +instance Random CFloat where + randomR = randomRFloating + random rng = case random rng of + (x,rng') -> (realToFrac (x::Float), rng') + +instance Random CDouble where + --randomR = randomRFloating + -- A MYSTERY: + -- Presently, this is showing better performance than the Double instance: + -- (And yet, if the Double instance uses randomFrac then its performance is much worse!) + --random = randomFrac + -- random rng = case random rng of + -- (x,rng') -> (realToFrac (x::Double), rng') + + + +randomBounded :: (RandomGen g, Random a, Bounded a) => g -> (a, g) +randomBounded = randomR (minBound, maxBound) + +-- The two integer functions below take an [inclusive,inclusive] range. +randomIvalIntegral :: (RandomGen g, Integral a) => (a, a) -> g -> (a, g) +randomIvalIntegral (l,h) = randomIvalInteger (toInteger l, toInteger h) + +-- {-# SPECIALIZE randomIvalInteger :: (Num a) => +-- (Integer, Integer) -> StdGen -> (a, StdGen) #-} + +randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g) +randomIvalInteger (l,h) rng + | l > h = randomIvalInteger (h,l) rng + | otherwise = case (f 1 0 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng') + where + (genlo, genhi) = genRange rng + b = fromIntegral genhi - fromIntegral genlo + 1 + + -- Probabilities of the most likely and least likely result + -- will differ at most by a factor of (1 +- 1/q). Assuming the RandomGen + -- is uniform, of course + + -- On average, log q / log b more random values will be generated + -- than the minimum + q = 1000 + k = h - l + 1 + magtgt = k * q + + -- generate random values until we exceed the target magnitude + f mag v g | mag >= magtgt = (v, g) + | otherwise = v' `seq`f (mag*b) v' g' where + (x,g') = next g + v' = (v * b + (fromIntegral x - fromIntegral genlo)) + + +{--- The continuous functions on the other hand take an [inclusive,exclusive) range. +randomFrac :: (RandomGen g, Fractional a) => g -> (a, g) +randomFrac = randomIvalDouble (0::Double,1) realToFrac-} + +--randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g) +--randomIvalDouble (l,h) fromDouble rng +-- | l > h = randomIvalDouble (h,l) fromDouble rng +-- | otherwise = +-- case (randomIvalInteger (toInteger (minBound::Int32), toInteger (maxBound::Int32)) rng) of +-- (x, rng') -> +-- let +-- scaled_x = +-- fromDouble (0.5*l + 0.5*h) + -- previously (l+h)/2, overflowed +-- fromDouble ((0.5*h - 0.5*l) / (0.5 * realToFrac int32Count)) * -- avoid overflow +-- fromIntegral (x::Int32) +-- in +-- (scaled_x, rng') + + +-- The global random number generator + +{- $globalrng #globalrng# + +There is a single, implicit, global random number generator of type +'StdGen', held in some global variable maintained by the 'IO' monad. It is +initialised automatically in some system-dependent fashion, for example, by +using the time of day, or Linux's kernel random number generator. To get +deterministic behaviour, use 'setStdGen'. +-} +{- +-- |Sets the global random number generator. +setStdGen :: StdGen -> IO () +setStdGen sgen = writeIORef theStdGen sgen + +-- |Gets the global random number generator. +getStdGen :: IO StdGen +getStdGen = readIORef theStdGen + +theStdGen :: IORef StdGen +theStdGen = unsafePerformIO $ do + rng <- mkStdRNG 0 + newIORef rng + +-- |Applies 'split' to the current global random generator, +-- updates it with one of the results, and returns the other. +newStdGen :: IO StdGen +newStdGen = atomicModifyIORef' theStdGen split + +{- |Uses the supplied function to get a value from the current global +random generator, and updates the global generator with the new generator +returned by the function. For example, @rollDice@ gets a random integer +between 1 and 6: + +> rollDice :: IO Int +> rollDice = getStdRandom (randomR (1,6)) + +-} + +getStdRandom :: (StdGen -> (a,StdGen)) -> IO a +getStdRandom f = atomicModifyIORef' theStdGen (swap . f) + where swap (v,g) = (g,v) +-} diff --git a/src/System/Random/PCG/Internal.hs b/src/System/Random/PCG32/Internal.hs similarity index 99% rename from src/System/Random/PCG/Internal.hs rename to src/System/Random/PCG32/Internal.hs index e828840f2..aae6ce0ed 100644 --- a/src/System/Random/PCG/Internal.hs +++ b/src/System/Random/PCG32/Internal.hs @@ -6,7 +6,7 @@ -- Standard PCG32 Random Number Generator with chosen streams, written in -- pure haskell. See for details. -- -module System.Random.PCG.Internal +module System.Random.PCG32.Internal ( -- * PCG 32 PCG32 (..) , seed diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index 14eb41b82..ef04cc4fe 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -1,11 +1,13 @@ {-# LANGUAGE ScopedTypeVariables, BangPatterns, UnboxedTuples, MagicHash, GADTs #-} -{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE DeriveFunctor, DeriveDataTypeable #-} module System.Random.SplitMix.Internal( nextSeedSplitMix ,splitGeneratorSplitMix ,nextWord64SplitMix + ,splitGeneratorSplitMix# + ,nextWord64SplitMix# ,SplitMix64(..) ,Random(..) ,sampleWord64Random @@ -17,6 +19,7 @@ import qualified Data.Bits as DB import Data.Bits (xor,(.|.)) import Data.Word(Word64) import Data.Functor.Identity +import Data.Data(Data(),Typeable()) {-# SPECIALIZE popCount :: Word64 -> Word64 #-} {-# SPECIALIZE popCount :: Int -> Word64 #-} @@ -73,6 +76,7 @@ type SplitMix64 = (# Word64# , Word64# #) data SplitMix64 = SplitMix64 { sm64seed :: {-# UNPACK #-} !Word64 ,sm64Gamma :: {-# UNPACK #-} !Word64 } + deriving (Eq,Read,Show,Data,Typeable) advanceSplitMix :: SplitMix64 -> SplitMix64 @@ -96,7 +100,7 @@ instance Functor Random where instance Applicative Random where pure = \ x -> Random# $ \ s -> (# x , s #) (<*>) = \ (Random# frmb) (Random# rma) -> Random# $ \ s -> - let (# fseed, maseed #) = splitGeneratorSplitMix s + let (# fseed, maseed #) = splitGeneratorSplitMix# s (# f , _boringSeed #) = frmb fseed (# a , newSeed #) = rma maseed in (# f a , newSeed #) @@ -105,12 +109,12 @@ instance Monad Random where (>>=) = \(Random# ma) f -> Random# $ \ s -> - let (# splitSeed , nextSeed #) = splitGeneratorSplitMix s + let (# splitSeed , nextSeed #) = splitGeneratorSplitMix# s (# a, _boringSeed #) = ma splitSeed in unRandom# (f a) nextSeed sampleWord64Random :: Random Word64 -sampleWord64Random = Random# nextWord64SplitMix +sampleWord64Random = Random# nextWord64SplitMix# newtype RandomT m a = RandomT# { unRandomT# :: (SplitMix64 -> m (a , SplitMix64) ) } @@ -122,7 +126,7 @@ instance Functor m => Functor (RandomT m) where instance Applicative m => Applicative (RandomT m) where pure = \ x -> RandomT# $ \ s -> pure ( x , s ) (<*>) = \ (RandomT# frmb) (RandomT# rma) -> RandomT# $ \ s -> - let (# !fseed, !maseed #) = splitGeneratorSplitMix s + let (# !fseed, !maseed #) = splitGeneratorSplitMix# s mfOldSeed= frmb fseed mArgNewSeed = rma maseed in (fmap (\(f,_s)-> \(x,newSeed)-> (f x, newSeed) ) mfOldSeed) @@ -131,7 +135,7 @@ instance Applicative m => Applicative (RandomT m) where instance Monad m => Monad (RandomT m) where (>>=) = \ (RandomT# ma) f -> RandomT# $ \ s -> - let (# fseed, nextSeed #) = splitGeneratorSplitMix s + let (# fseed, nextSeed #) = splitGeneratorSplitMix# s in do (a,_boring) <- ma fseed @@ -139,25 +143,36 @@ instance Monad m => Monad (RandomT m) where sampleWord64RandomT :: Applicative m => RandomT m Word64 sampleWord64RandomT = RandomT# $ \ s -> - let (# !w, !ngen #) = nextWord64SplitMix s + let (# !w, !ngen #) = nextWord64SplitMix# s in pure (w, ngen) --instance PrimMonad m => PrimMonad (RandomT m) where -- primitive = \ m -> -- {-# INLINE #-} -nextWord64SplitMix :: SplitMix64 -> (# Word64 , SplitMix64 #) -nextWord64SplitMix gen = mixedRes `seq` (# mixedRes , newgen #) +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 #)) +{-# INLINE nextWord64SplitMix #-} +nextWord64SplitMix :: SplitMix64 -> ( Word64 , SplitMix64 ) +nextWord64SplitMix gen = (mixedRes, newgen) where - (# splitSeed , nextGen #) = nextWord64SplitMix gen + (# mixedRes,newgen #) = nextWord64SplitMix# 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 - +{-# INLINE splitGeneratorSplitMix #-} +splitGeneratorSplitMix :: SplitMix64 -> (SplitMix64 , SplitMix64) +splitGeneratorSplitMix gen = (splitGen, nextNextGen) + where + (# splitGen, nextNextGen #) = splitGeneratorSplitMix# gen diff --git a/src/System/Random/SplitMix/Internal/Splitting.hs b/src/System/Random/SplitMix/Internal/Splitting.hs new file mode 100644 index 000000000..be6b5669f --- /dev/null +++ b/src/System/Random/SplitMix/Internal/Splitting.hs @@ -0,0 +1,2 @@ +module System.Random.SplitMix.Internal.Splitting where + From c6c3ed72480bf35c90dfdecf8d03ea9d5da64e7c Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Mon, 25 Dec 2017 20:13:33 -0500 Subject: [PATCH 29/31] Notes / References for constructing a splitmix rng that follows the lea&steele fixes that weren't in the paper but are in the subsequent JDK releases, and are also noted by the PCG algorithm author Melissa ONeil --- src/System/Random/SplitMix/Internal.hs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index ef04cc4fe..c40e45196 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -21,6 +21,20 @@ import Data.Word(Word64) import Data.Functor.Identity import Data.Data(Data(),Typeable()) +{- +splitmix constants follows +https://github.com/dmlloyd/openjdk/blob/67672eec97164de10a9ca83ddbcef6b42816ed04/src/java.base/share/classes/java/util/SplittableRandom.java + +see also +http://hg.openjdk.java.net/jdk/jdk10/file/bffcbf07ea88/src/java.base/share/classes/java/util/SplittableRandom.java + +ie the variant found in JDK >=8 + +see also discussion on the melissa o'neil pcg blog about +splitmix +http://www.pcg-random.org/posts/bugs-in-splitmix.html +-} + {-# SPECIALIZE popCount :: Word64 -> Word64 #-} {-# SPECIALIZE popCount :: Int -> Word64 #-} {-# SPECIALIZE popCount :: Word -> Word64 #-} From 22a2a16bd62edd553b4f7f2e9eedc26cbf8850d8 Mon Sep 17 00:00:00 2001 From: Carter Tazio Schonwald Date: Thu, 1 Mar 2018 13:51:20 -0500 Subject: [PATCH 30/31] add ord insance to splitmix generator datatype --- src/System/Random/SplitMix/Internal.hs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs index c40e45196..b6b9a79e6 100644 --- a/src/System/Random/SplitMix/Internal.hs +++ b/src/System/Random/SplitMix/Internal.hs @@ -90,7 +90,7 @@ type SplitMix64 = (# Word64# , Word64# #) data SplitMix64 = SplitMix64 { sm64seed :: {-# UNPACK #-} !Word64 ,sm64Gamma :: {-# UNPACK #-} !Word64 } - deriving (Eq,Read,Show,Data,Typeable) + deriving (Eq,Ord,Read,Show,Data,Typeable) advanceSplitMix :: SplitMix64 -> SplitMix64 From 78f70ad25b705371b4a2c82255624dc2fc0d97a7 Mon Sep 17 00:00:00 2001 From: Ron Wolf Date: Sat, 1 Sep 2018 14:16:34 -0400 Subject: [PATCH 31/31] Insert build status image inline --- readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/readme.md b/readme.md index e87917c42..61de059f6 100644 --- a/readme.md +++ b/readme.md @@ -1,4 +1,4 @@ -https://travis-ci.org/cartazio/random.svg?branch=master +![Build status](https://travis-ci.org/cartazio/random.svg?branch=master) # Random