From 326cc7edbc318c25d3f38d424ff7ff4bf1617475 Mon Sep 17 00:00:00 2001 From: Leonhard Markert Date: Wed, 13 May 2020 11:00:10 +0200 Subject: [PATCH] master -> v1.1 --- .darcs-boring | 5 + .gitignore | 13 +- .travis.yml | 108 +----- Benchmark/BinSearch.hs | 142 ++++++++ Benchmark/Makefile | 24 ++ Benchmark/SimpleRNGBench.hs | 322 +++++++++++++++++ CHANGELOG.md | 27 +- DEVLOG.md | 196 +++++++++++ LICENSE | 83 +++-- README.md | 18 + Setup.hs | 4 + {src/System => System}/Random.hs | 328 +++++++++++------- TODO | 15 - cmmbits/floatsAndBits.cmm | 68 ---- prologue.txt | 1 + random.cabal | 168 ++++----- readme.md | 5 - src/Data/Distribution/FloatingInterval.hs | 101 ------ src/Data/Distribution/Integers.hs | 50 --- src/Data/Distribution/Normal.hs | 36 -- src/Data/Distribution/Permutation.hs | 72 ---- src/Data/Random/Utils.hs | 136 -------- src/System/Random/PCG32/Internal.hs | 160 --------- src/System/Random/SplitMix/Internal.hs | 192 ---------- .../Random/SplitMix/Internal/Splitting.hs | 2 - testCast/WordFloat.hs | 31 -- tests/Makefile | 14 + tests/T7936.hs | 14 + tests/TestRandomIOs.hs | 20 ++ tests/TestRandomRs.hs | 22 ++ tests/all.T | 10 + tests/random1283.hs | 40 +++ tests/rangeTest.hs | 135 +++++++ tests/rangeTest.stdout | 67 ++++ tests/slowness.hs | 55 +++ 35 files changed, 1460 insertions(+), 1224 deletions(-) create mode 100644 .darcs-boring create mode 100644 Benchmark/BinSearch.hs create mode 100644 Benchmark/Makefile create mode 100644 Benchmark/SimpleRNGBench.hs create mode 100644 DEVLOG.md create mode 100644 README.md rename {src/System => System}/Random.hs (69%) delete mode 100644 TODO delete mode 100644 cmmbits/floatsAndBits.cmm create mode 100644 prologue.txt delete mode 100644 readme.md delete mode 100644 src/Data/Distribution/FloatingInterval.hs delete mode 100644 src/Data/Distribution/Integers.hs delete mode 100644 src/Data/Distribution/Normal.hs delete mode 100644 src/Data/Distribution/Permutation.hs delete mode 100644 src/Data/Random/Utils.hs delete mode 100644 src/System/Random/PCG32/Internal.hs delete mode 100644 src/System/Random/SplitMix/Internal.hs delete mode 100644 src/System/Random/SplitMix/Internal/Splitting.hs delete mode 100644 testCast/WordFloat.hs create mode 100644 tests/Makefile create mode 100644 tests/T7936.hs create mode 100644 tests/TestRandomIOs.hs create mode 100644 tests/TestRandomRs.hs create mode 100644 tests/all.T create mode 100644 tests/random1283.hs create mode 100644 tests/rangeTest.hs create mode 100644 tests/rangeTest.stdout create mode 100644 tests/slowness.hs diff --git a/.darcs-boring b/.darcs-boring new file mode 100644 index 000000000..6682606a0 --- /dev/null +++ b/.darcs-boring @@ -0,0 +1,5 @@ +^dist(/|$) +^setup(/|$) +^GNUmakefile$ +^Makefile.local$ +^.depend(.bak)?$ diff --git a/.gitignore b/.gitignore index c70beef85..41c6d8c61 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,12 @@ -cabal.project.local +*~ + +Thumbs.db +.DS_Store + +GNUmakefile +dist-install/ +ghc.mk + +dist +.cabal-sandbox +cabal.sandbox.config diff --git a/.travis.yml b/.travis.yml index d3e6cac1d..6a03bcb12 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,103 +1,5 @@ -# 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 +language: haskell +ghc: + - 7.4 + - 7.6 + - 7.8 diff --git a/Benchmark/BinSearch.hs b/Benchmark/BinSearch.hs new file mode 100644 index 000000000..f61164855 --- /dev/null +++ b/Benchmark/BinSearch.hs @@ -0,0 +1,142 @@ + +{- + 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 new file mode 100644 index 000000000..8a84e6479 --- /dev/null +++ b/Benchmark/Makefile @@ -0,0 +1,24 @@ + + +#OPTS= -O2 -Wall -XCPP +OPTS= -O2 -Wall -XCPP -Werror + +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/SimpleRNGBench.hs b/Benchmark/SimpleRNGBench.hs new file mode 100644 index 000000000..c25b75d47 --- /dev/null +++ b/Benchmark/SimpleRNGBench.hs @@ -0,0 +1,322 @@ +{-# 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.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 :: Integral a => a -> String +commaint n = + reverse $ concat $ + intersperse "," $ + chunk 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 + +-- gen_mwc <- create + 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 71040582c..15c882af9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,26 @@ -# Revision history for random +# 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 -## 2.0.0.0 -- YYYY-mm-dd +# 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 -* First version. Released on an unsuspecting world. diff --git a/DEVLOG.md b/DEVLOG.md new file mode 100644 index 000000000..6e0b28dc4 --- /dev/null +++ b/DEVLOG.md @@ -0,0 +1,196 @@ + + +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 0c53b4a83..06bb64148 100644 --- a/LICENSE +++ b/LICENSE @@ -1,26 +1,63 @@ -Copyright (c) 2017, Carter Tazio Schonwald +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. All rights reserved. Redistribution and use in source and binary forms, with or without -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. +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. + +----------------------------------------------------------------------------- diff --git a/README.md b/README.md new file mode 100644 index 000000000..9d5bb51b2 --- /dev/null +++ b/README.md @@ -0,0 +1,18 @@ +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] +(http://www.haskell.org/ghc/docs/latest/html/libraries/haskell98/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 9a994af67..6fa548caf 100644 --- a/Setup.hs +++ b/Setup.hs @@ -1,2 +1,6 @@ +module Main (main) where + import Distribution.Simple + +main :: IO () main = defaultMain diff --git a/src/System/Random.hs b/System/Random.hs similarity index 69% rename from src/System/Random.hs rename to System/Random.hs index f6bde7f35..ab7727405 100644 --- a/src/System/Random.hs +++ b/System/Random.hs @@ -1,4 +1,3 @@ -{-# LANGUAGE CPP #-} #if __GLASGOW_HASKELL__ >= 701 {-# LANGUAGE Trustworthy #-} #endif @@ -8,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 @@ -19,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. @@ -41,38 +40,40 @@ #include "MachDeps.h" module System.Random - ( - - -- $intro - - -- * Random number generators + ( + -- $intro - RandomGen(next, genRange) - , SplittableGen(split) + -- * Random number generators - -- ** Standard random number generators - --, StdGen - --, mkStdGen +#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 + -- ** 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 @@ -81,19 +82,24 @@ 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 ) +#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' ) +import Data.IORef ( atomicModifyIORef' ) #else import Data.IORef ( atomicModifyIORef ) #endif ---import Numeric ( readDec ) +import Numeric ( readDec ) #ifdef __GLASGOW_HASKELL__ import GHC.Exts ( build ) @@ -113,19 +119,36 @@ atomicModifyIORef' ref f = do 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 -> (Word64, g) + next :: g -> (Int, g) -- |The 'genRange' operation yields the range of values returned by -- the generator. @@ -144,16 +167,16 @@ class RandomGen g where -- a different range to the generator passed to 'next'. -- -- The default definition spans the full range of 'Int'. - genRange :: g -> (Word64,Word64) + 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 @@ -178,30 +201,31 @@ 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 --- = StdGen !Int32 !Int32 -{- +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 . + showsPrec p (StdGen s1 s2) = + showsPrec p s1 . showChar ' ' . showsPrec p s2 @@ -210,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 @@ -222,15 +246,15 @@ 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) - | +{- | 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 @@ -242,15 +266,14 @@ 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 --} {- | With a source of random number supply in hand, the 'Random' class allows the @@ -295,12 +318,12 @@ class Random a where -- | 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) + 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 + randomIO = getStdRandom random -- | Produce an infinite list-equivalent of random values. {-# INLINE buildRandoms #-} @@ -317,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 @@ -325,9 +348,10 @@ instance Random Int16 where randomR = randomIvalIntegral; random = randomBo 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 @@ -354,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 = - -- case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of - -- (x,g') -> (chr x, g') - --random g = randomR (minBound,maxBound) 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 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 @@ -368,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 @@ -412,19 +436,23 @@ 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 + 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') - + 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) @@ -433,9 +461,9 @@ randomBounded = randomR (minBound, maxBound) 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) #-} - +{-# 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 @@ -454,31 +482,66 @@ 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 v' = (v * b + (fromIntegral x - fromIntegral genlo)) -{--- The continuous functions on the other hand take an [inclusive,exclusive) range. +-- 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') - +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 @@ -490,7 +553,7 @@ 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 @@ -522,4 +585,25 @@ between 1 and 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/TODO b/TODO deleted file mode 100644 index ef0dea00b..000000000 --- a/TODO +++ /dev/null @@ -1,15 +0,0 @@ -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/cmmbits/floatsAndBits.cmm b/cmmbits/floatsAndBits.cmm deleted file mode 100644 index acfa456be..000000000 --- a/cmmbits/floatsAndBits.cmm +++ /dev/null @@ -1,68 +0,0 @@ -#include "Cmm.h" -#include "MachDeps.h" - -#if WORD_SIZE_IN_BITS == 64 -#define DOUBLE_SIZE_WDS 1 -#else -#define DOUBLE_SIZE_WDS 2 -#endif - -cts_random_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); -} - -cts_random_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); -} - -cts_random_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); -} - -cts_random_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/prologue.txt b/prologue.txt new file mode 100644 index 000000000..ea0344b9c --- /dev/null +++ b/prologue.txt @@ -0,0 +1 @@ +Random number library. diff --git a/random.cabal b/random.cabal index 2d4f18d21..fd29840fb 100644 --- a/random.cabal +++ b/random.cabal @@ -1,110 +1,70 @@ --- Initial random.cabal generated by cabal init. For further --- documentation, see http://haskell.org/cabal/users-guide/ +name: random +version: 1.1 --- 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: 1.3.0.0 --- A short (one-line) description of the package. -synopsis: Random number generation library for haskell --- A longer description of the package. --- description: +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. --- 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 - .travis.yml - readme.md - - - --- Constraint on the version of Cabal needed to build this package. -cabal-version: >=1.10 - - -library - -- Modules exported` by the library. - exposed-modules: - System.Random - System.Random.SplitMix.Internal - System.Random.SplitMix.Internal.Splitting - System.Random.PCG32.Internal - Data.Distribution.FloatingInterval - Data.Distribution.Normal - Data.Distribution.Permutation - Data.Distribution.Integers - Data.Random.Utils - - c-sources: cmmbits/floatsAndBits.cmm - -- 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.8 && <4.12 - ,ghc-prim - ,entropy == 0.3.* - ,numeric-extras == 0.1.* - ,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. - hs-source-dirs: src - - ghc-options: -O2 -Wall - - -- 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 + .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: http://git.haskell.org/packages/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 diff --git a/readme.md b/readme.md deleted file mode 100644 index e87917c42..000000000 --- a/readme.md +++ /dev/null @@ -1,5 +0,0 @@ -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 deleted file mode 100644 index 62cda7514..000000000 --- a/src/Data/Distribution/FloatingInterval.hs +++ /dev/null @@ -1,101 +0,0 @@ - -{- | 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/Integers.hs b/src/Data/Distribution/Integers.hs deleted file mode 100644 index e0d4959e5..000000000 --- a/src/Data/Distribution/Integers.hs +++ /dev/null @@ -1,50 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -module Data.Distribution.Integers where - -import Data.Word(Word64) -import Data.Bits - -{- -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 --} - -{- | @'sampleWordRange' wordSampler (lo,hi)@ will return a uniform sample from the closed interval -@[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 - | 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 - -- 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)) - - diff --git a/src/Data/Distribution/Normal.hs b/src/Data/Distribution/Normal.hs deleted file mode 100644 index 910db8cdc..000000000 --- a/src/Data/Distribution/Normal.hs +++ /dev/null @@ -1,36 +0,0 @@ -{-# 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) - - diff --git a/src/Data/Distribution/Permutation.hs b/src/Data/Distribution/Permutation.hs deleted file mode 100644 index 9bd642638..000000000 --- a/src/Data/Distribution/Permutation.hs +++ /dev/null @@ -1,72 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables, GADTs#-} -module Data.Distribution.Permutation where - -import Control.Monad.Primitive as Prim ---import Data.Primitive.Array as DPA -import Data.Word(Word32) -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. - -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 - -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 - 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 :: 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)] <- 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)] - -> 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" - | otherwise = do swapSpots 0 ls ; DV.unsafeFreeze marr - where - arrayLength :: Int - 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 - 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 DVM.swap marr from to - swapSpots (ix +1) rest - - diff --git a/src/Data/Random/Utils.hs b/src/Data/Random/Utils.hs deleted file mode 100644 index 4bd2ec1ad..000000000 --- a/src/Data/Random/Utils.hs +++ /dev/null @@ -1,136 +0,0 @@ -{-# LANGUAGE Trustworthy #-} -{-# LANGUAGE CPP - , GHCForeignImportPrim - , MagicHash - , UnboxedTuples - , UnliftedFFITypes - #-} -{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} - -#include "MachDeps.h" -module Data.Random.Utils( - castWord64ToDouble - ,castDoubleToWord64 - ,castWord32ToFloat - ,castFloatToWord32 - ,RandomCastIEEE(..)) 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 -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 RandomCastIEEE word ieee | word -> ieee , ieee -> word where - toIEEE :: word -> ieee - fromIEEE :: ieee -> word - -instance RandomCastIEEE Word32 Float where - {-# INLINE toIEEE #-} - {-# INLINE fromIEEE #-} - toIEEE = castWord32ToFloat - fromIEEE = castFloatToWord32 -instance RandomCastIEEE 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 "cts_random_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 "cts_random_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 "cts_random_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 "cts_random_stg_doubleToWord64zhPrivate" -#if WORD_SIZE_IN_BITS == 64 - stgDoubleToWord64 :: Double# -> Word# -#else - stgDoubleToWord64 :: Double# -> Word64# -#endif diff --git a/src/System/Random/PCG32/Internal.hs b/src/System/Random/PCG32/Internal.hs deleted file mode 100644 index aae6ce0ed..000000000 --- a/src/System/Random/PCG32/Internal.hs +++ /dev/null @@ -1,160 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE DeriveDataTypeable #-} -{-# LANGUAGE DeriveGeneric #-} --- | --- --- Standard PCG32 Random Number Generator with chosen streams, written in --- pure haskell. See for details. --- -module System.Random.PCG32.Internal - ( -- * PCG 32 - PCG32 (..) - , seed - , newPCG32 - - -- * Generating random numbers - , next32 - , next64 - , bounded - - -- * Generator utilities - , advancePCG32 - , split - ) where - -import Data.Bits -import Data.Data -import Foreign -import GHC.Generics - --- | The multiple sequence varient of the pcg random number generator. -data PCG32 = PCG32 - {-# UNPACK #-} !Word64 -- state - {-# UNPACK #-} !Word64 -- inc - deriving (Show, Ord, Eq, Data, Typeable, Generic) - --- | Create a new generator from two words. -newPCG32 :: Word64 -> Word64 -> PCG32 -newPCG32 = \a b -> - let s = state (PCG32 (a + i) i) - i = (b `shiftL` 1) .|. 1 - in PCG32 s i -{-# INLINE newPCG32 #-} - --- | Fixed seed. -seed :: PCG32 -seed = PCG32 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 :: PCG32 -> Word64 -state (PCG32 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 :: 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 :: PCG32 -> (Word32, PCG32) -next32 = \g@(PCG32 _ inc) -> - let Pair s' r = pair g - in (r, PCG32 s' inc) -{-# INLINE next32 #-} - --- | Return a random 'Word64' with the generator incremented by two steps. -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 -> 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 (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 -> PCG32 -> (Word32, PCG32) -bounded = \b g@(PCG32 _ inc) -> - let !(Pair s1 w) = bounded' b g - in (w, PCG32 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 :: 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 = 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 - :: 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. -advancePCG32 :: Word64 -> PCG32 -> PCG32 -advancePCG32 = \d (PCG32 s inc) -> PCG32 (advancing d s multiplier inc) inc -{-# INLINE advancePCG32 #-} - --- Misc ---------------------------------------------------------------- - -wordsTo64Bit :: Integral a => Word32 -> Word32 -> a -wordsTo64Bit x y = - fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64) -{-# INLINE wordsTo64Bit #-} - diff --git a/src/System/Random/SplitMix/Internal.hs b/src/System/Random/SplitMix/Internal.hs deleted file mode 100644 index b6b9a79e6..000000000 --- a/src/System/Random/SplitMix/Internal.hs +++ /dev/null @@ -1,192 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables, BangPatterns, UnboxedTuples, MagicHash, GADTs #-} -{-# LANGUAGE DeriveFunctor, DeriveDataTypeable #-} - - -module System.Random.SplitMix.Internal( - nextSeedSplitMix - ,splitGeneratorSplitMix - ,nextWord64SplitMix - ,splitGeneratorSplitMix# - ,nextWord64SplitMix# - ,SplitMix64(..) - ,Random(..) - ,sampleWord64Random - ,RandomT(..) - ,sampleWord64RandomT - ) where - -import qualified Data.Bits as DB -import Data.Bits (xor,(.|.)) -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 #-} -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 - -{- - -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 } - deriving (Eq,Ord,Read,Show,Data,Typeable) - - -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 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 #) - -instance Monad Random where - (>>=) = - \(Random# ma) f -> - Random# $ \ s -> - let (# splitSeed , nextSeed #) = splitGeneratorSplitMix# s - (# 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) ) } - -instance Functor m => Functor (RandomT m) where - fmap = \ f (RandomT# mf) -> - RandomT# $ \ seed -> - fmap (\(a,s) -> (f a, s) ) $ mf seed - -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 - 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) f -> RandomT# $ \ s -> - let (# fseed, nextSeed #) = splitGeneratorSplitMix# s - in - do - (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) - ---instance PrimMonad m => PrimMonad (RandomT m) where --- primitive = \ m -> --- {-# INLINE #-} - -nextWord64SplitMix# :: SplitMix64 -> (# Word64 , SplitMix64 #) -nextWord64SplitMix# gen = mixedRes `seq` (# mixedRes , newgen #) - where - mixedRes = mix64 premixres - (# premixres , newgen #) = nextSeedSplitMix gen - -{-# INLINE nextWord64SplitMix #-} -nextWord64SplitMix :: SplitMix64 -> ( Word64 , SplitMix64 ) -nextWord64SplitMix gen = (mixedRes, newgen) - where - (# 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 deleted file mode 100644 index be6b5669f..000000000 --- a/src/System/Random/SplitMix/Internal/Splitting.hs +++ /dev/null @@ -1,2 +0,0 @@ -module System.Random.SplitMix.Internal.Splitting where - diff --git a/testCast/WordFloat.hs b/testCast/WordFloat.hs deleted file mode 100644 index b4dea2179..000000000 --- a/testCast/WordFloat.hs +++ /dev/null @@ -1,31 +0,0 @@ - -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 )) - ]] - - diff --git a/tests/Makefile b/tests/Makefile new file mode 100644 index 000000000..39c71491b --- /dev/null +++ b/tests/Makefile @@ -0,0 +1,14 @@ +# 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 new file mode 100644 index 000000000..cfea911bc --- /dev/null +++ b/tests/T7936.hs @@ -0,0 +1,14 @@ +-- 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 new file mode 100644 index 000000000..d8a00cc7c --- /dev/null +++ b/tests/TestRandomIOs.hs @@ -0,0 +1,20 @@ +-- 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 new file mode 100644 index 000000000..cdae106a6 --- /dev/null +++ b/tests/TestRandomRs.hs @@ -0,0 +1,22 @@ +-- 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 new file mode 100644 index 000000000..f1675ed5c --- /dev/null +++ b/tests/all.T @@ -0,0 +1,10 @@ + +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 new file mode 100644 index 000000000..33fc48862 --- /dev/null +++ b/tests/random1283.hs @@ -0,0 +1,40 @@ +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 new file mode 100644 index 000000000..ac62c71dd --- /dev/null +++ b/tests/rangeTest.hs @@ -0,0 +1,135 @@ +{-# 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 new file mode 100644 index 000000000..55ccaffb4 --- /dev/null +++ b/tests/rangeTest.stdout @@ -0,0 +1,67 @@ +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 new file mode 100644 index 000000000..162a4cff0 --- /dev/null +++ b/tests/slowness.hs @@ -0,0 +1,55 @@ + +-- 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.