mwc-random 0.14.0.0 → 0.15.0.0

raw patch · 8 files changed

+618/−182 lines, 8 filesdep +QuickCheckdep +bytestringdep +doctestdep ~basedep ~primitivedep ~vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, bytestring, doctest, gauge, mersenne-random, mwc-random, random, tasty, tasty-hunit, tasty-quickcheck

Dependency ranges changed: base, primitive, vector

API changes (from Hackage documentation)

Files

• System/Random/MWC.hs +254/−126
• System/Random/MWC/CondensedTable.hs +3/−6
• System/Random/MWC/Distributions.hs +45/−47
• bench/Benchmark.hs +111/−0
• changelog.md +19/−0
• mwc-random.cabal +62/−3
• tests/doctests.hs +4/−0
• tests/props.hs +120/−0

@@ -1,5 +1,6 @@ {-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, FlexibleContexts,-    MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples+    FlexibleInstances, MultiParamTypeClasses, MagicHash, Rank2Types,+    ScopedTypeVariables, TypeFamilies, UnboxedTuples     #-} -- | -- Module    : System.Random.MWC@@ -10,68 +11,123 @@ -- Stability   : experimental -- Portability : portable ----- Pseudo-random number generation.  This module contains code for--- generating high quality random numbers that follow a uniform--- distribution.+-- Pseudo-random number generation using Marsaglia's MWC256, (also+-- known as MWC8222) multiply-with-carry generator, which has a period+-- of \(2^{8222}\) and fares well in tests of randomness.  It is also+-- extremely fast, between 2 and 3 times faster than the Mersenne+-- Twister. There're two representation of generator: 'Gen' which is+-- generator that uses in-place mutation and 'Seed' which is immutable+-- snapshot of generator's state. ----- For non-uniform distributions, see the--- 'System.Random.MWC.Distributions' module. ----- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222)--- multiply-with-carry generator, which has a period of 2^8222 and--- fares well in tests of randomness.  It is also extremely fast,--- between 2 and 3 times faster than the Mersenne Twister.+-- == Initialization ----- The generator state is stored in the 'Gen' data type. It can be--- created in several ways:+-- Generator could be initialized in several ways. One is to obtain+-- randomness from operating system using 'createSystemRandom',+-- 'createSystemSeed' or 'withSystemRandomST' (All examples assume+-- that @System.Random.Stateful@ is imported) -----   1. Using the 'withSystemRandom' call, which creates a random state.+-- >>> g <- createSystemRandom+-- >>> uniformM g :: IO Int+-- ... -----   2. Supply your own seed to 'initialize' function.+-- >>> withSystemRandomST $ \g -> uniformM g :: IO Int+-- ... -----   3. Finally, 'create' makes a generator from a fixed seed.---      Generators created in this way aren't really random.+-- Deterministically create generator from given seed using+-- 'initialize' function: ----- For repeatability, the state of the generator can be snapshotted--- and replayed using the 'save' and 'restore' functions.+-- >>> import qualified Data.Vector.Unboxed as U+-- >>> import System.Random.Stateful+-- >>> g <- initialize $ U.fromList [1,2,3]+-- >>> uniformRM (1,200) g :: IO Int+-- 101 ----- The simplest use is to generate a vector of uniformly distributed values:+-- Last way is to create generator with fixed seed which could be+-- useful in testing ----- @---   vs \<- 'withSystemRandom' . 'asGenST' $ \\gen -> 'uniformVector' gen 100--- @+-- >>> g <- create+-- >>> uniformM g :: IO Int+-- -8765701622605876598 ----- These values can be of any type which is an instance of the class--- 'Variate'. ----- To generate random values on demand, first 'create' a random number--- generator.+-- == Generation of random numbers ----- @---   gen <- 'create'--- @+-- Recommended way of generating random numbers in simple cases like+-- generating uniformly distributed random number in range or value+-- uniformly distributed in complete type domain is to use+-- 'UniformRange' and 'Uniform' type classes. Note that while small+-- self-contained examples usually require explicit annotations+-- usually result type could be inferred. ----- Hold onto this generator and use it wherever random values are--- required (creating a new generator is expensive compared to--- generating a random number, so you don't want to throw them--- away). Get a random value using 'uniform' or 'uniformR':+-- This example simulates 20 throws of fair 6-sided dice: ----- @---   v <- 'uniform' gen--- @+-- >>> g <- create+-- >>> replicateM 20 $ uniformRM (1, 6::Integer) g+-- [3,4,3,1,4,6,1,6,1,4,2,2,3,2,4,2,5,1,3,5] ----- @---   v <- 'uniformR' (1, 52) gen--- @+-- For generating full range of possible values one could use+-- 'uniformM'. This example generates 10 random bytes, or equivalently+-- 10 throws of 256-sided dice:+--+-- >>> g <- create+-- >>> replicateM 10 $ uniformM g :: IO [Word8]+-- [209,138,126,150,165,15,69,203,155,146]+--+-- There're special functions for generation of @Doubles@ and @Float+-- in unit interval: 'Random.uniformDouble01M',+-- 'Random.uniformDoublePositive01M', 'Random.uniformFloat01M',+-- 'Random.uniformFloatPositive01M':+--+-- >>> uniformDouble01M =<< create+-- 0.5248103628705498+-- >>> uniformFloat01M =<< create+-- 0.5248104+--+-- For normal distribution and others see modules+-- "System.Random.MWC.Distributions" and+-- "System.Random.MWC.CondensedTable". Note that they could be used+-- with any other generator implementing 'Random.StatefulGen' API+--+-- There're special cases for generating random vectors and+-- bytestrings. For example in order to generate random 10-byte+-- sequences as unboxed vector or bytestring:+--+-- >>> g <- create+-- >>> uniformVector g 10 :: IO (U.Vector Word8)+-- [209,138,126,150,165,15,69,203,155,146]+--+-- >>> import qualified Data.ByteString as BS+-- >>> g <- create+-- >>> BS.unpack <$> uniformByteStringM 10 g+-- [138,242,130,33,209,248,89,134,150,180]+--+-- Note that 'Random.uniformByteStringM' produces different result+-- from 'uniformVector' since it uses PRNG's output more efficently.+--+--+-- == State handling+--+-- For repeatability, the state of the generator can be snapshotted+-- and replayed using the 'save' and 'restore' functions. Following+-- example shows how to save and restore generator:+--+-- >>> g <- create+-- >>> replicateM_ 10 (uniformM g :: IO Word64)+-- >>> s <- save g+-- >>> uniformM g :: IO Word32+-- 1771812561+-- >>> uniformM =<< restore s :: IO Word32+-- 1771812561 module System.Random.MWC     (     -- * Gen: Pseudo-Random Number Generators       Gen     , create     , initialize-    , withSystemRandom+    , createSystemSeed     , createSystemRandom-+    , withSystemRandomST     -- ** Type helpers     -- $typehelp     , GenIO@@ -80,6 +136,8 @@     , asGenST      -- * Variates: uniformly distributed values+    , Random.Uniform(..)+    , Random.UniformRange(..)     , Variate(..)     , uniformVector @@ -89,38 +147,35 @@     , toSeed     , save     , restore-+    -- * Deprecated+    , withSystemRandom     -- * References     -- $references     ) where -#if defined(__GLASGOW_HASKELL__) && !defined(__HADDOCK__)-#include "MachDeps.h"-#endif- import Control.Monad           (ap, liftM, unless)-import Control.Monad.Primitive (PrimMonad, PrimBase, PrimState, unsafePrimToIO)-import Control.Monad.ST        (ST)+import Control.Monad.Primitive (PrimMonad, PrimBase, PrimState, unsafePrimToIO, unsafeSTToPrim)+import Control.Monad.ST        (ST,runST) import Data.Bits               ((.&.), (.|.), shiftL, shiftR, xor) import Data.Int                (Int8, Int16, Int32, Int64)-import Data.IORef              (atomicModifyIORef, newIORef)+import Data.IORef              (IORef, atomicModifyIORef, newIORef) import Data.Typeable           (Typeable) import Data.Vector.Generic     (Vector) import Data.Word import qualified Data.Vector.Generic         as G+import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Unboxed         as I import qualified Data.Vector.Unboxed.Mutable as M import System.IO        (hPutStrLn, stderr) import System.IO.Unsafe (unsafePerformIO) import qualified Control.Exception as E-#if defined(mingw32_HOST_OS)-import Foreign.Ptr-import Foreign.C.Types-#endif import System.Random.MWC.SeedSource-+import qualified System.Random.Stateful as Random --- | The class of types for which we can generate uniformly+-- | NOTE: Consider use of more principled type classes+-- 'Random.Uniform' and 'Random.UniformRange' instead.+--+-- The class of types for which we can generate uniformly -- distributed random variates. -- -- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222)@@ -244,14 +299,6 @@     {-# INLINE uniform  #-}     {-# INLINE uniformR #-} -{--instance Variate Integer where-    uniform g = do-      u <- uniform g-      return $! fromIntegral (u :: Int)-    {-# INLINE uniform #-}--}- instance (Variate a, Variate b) => Variate (a,b) where     uniform g = (,) `liftM` uniform g `ap` uniform g     uniformR ((x1,y1),(x2,y2)) g = (,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g@@ -385,12 +432,36 @@ {-# INLINE initialize #-}  -- | An immutable snapshot of the state of a 'Gen'.-newtype Seed = Seed {-    -- | Convert seed into vector.-    fromSeed :: I.Vector Word32-    }-    deriving (Eq, Show, Typeable)+newtype Seed = Seed (I.Vector Word32)+  deriving (Eq, Show, Typeable) +-- | Convert seed into vector.+fromSeed :: Seed -> I.Vector Word32+fromSeed (Seed v) = v++-- | @since 0.15.0.0+instance (s ~ PrimState m, PrimMonad m) => Random.StatefulGen (Gen s) m where+  uniformWord32R u = uniformR (0, u)+  {-# INLINE uniformWord32R #-}+  uniformWord64R u = uniformR (0, u)+  {-# INLINE uniformWord64R #-}+  uniformWord8 = uniform+  {-# INLINE uniformWord8 #-}+  uniformWord16 = uniform+  {-# INLINE uniformWord16 #-}+  uniformWord32 = uniform+  {-# INLINE uniformWord32 #-}+  uniformWord64 = uniform+  {-# INLINE uniformWord64 #-}+  uniformShortByteString n g = unsafeSTToPrim (Random.genShortByteStringST n (uniform g))+  {-# INLINE uniformShortByteString #-}++-- | @since 0.15.0.0+instance PrimMonad m => Random.FrozenGen Seed m where+  type MutableGen Seed m = Gen (PrimState m)+  thawGen = restore+  freezeGen = save+ -- | Convert vector to 'Seed'. It acts similarily to 'initialize' and -- will accept any vector. If you want to pass seed immediately to -- restore you better call initialize directly since following law holds:@@ -410,33 +481,74 @@ {-# INLINE restore #-}  --- | Seed a PRNG with data from the system's fast source of--- pseudo-random numbers (\"@\/dev\/urandom@\" on Unix-like systems or--- @RtlGenRandom@ on Windows), then run the given action.+-- $seeding ----- This is a somewhat expensive function, and is intended to be called--- only occasionally (e.g. once per thread).  You should use the `Gen`--- it creates to generate many random numbers.-withSystemRandom :: PrimBase m-                 => (Gen (PrimState m) -> m a) -> IO a-withSystemRandom act = do-  seed <- acquireSeedSystem 256 `E.catch` \(_::E.IOException) -> do-    seen <- atomicModifyIORef warned ((,) True)+-- Library provides several functions allowing to intialize generator+-- using OS-provided randomness: \"@\/dev\/urandom@\" on Unix-like+-- systems or @RtlGenRandom@ on Windows. This is a somewhat expensive+-- function, and is intended to be called only occasionally (e.g. once+-- per thread).  You should use the `Gen` it creates to generate many+-- random numbers.++createSystemRandomList :: IO [Word32]+createSystemRandomList = do+  acquireSeedSystem 256 `E.catch` \(_::E.IOException) -> do+    seen <- atomicModifyIORef seedCreatetionWarned ((,) True)     unless seen $ E.handle (\(_::E.IOException) -> return ()) $ do       hPutStrLn stderr $ "Warning: Couldn't use randomness source " ++ randomSourceName       hPutStrLn stderr ("Warning: using system clock for seed instead " ++                         "(quality will be lower)")     acquireSeedTime-  unsafePrimToIO $ initialize (I.fromList seed) >>= act-  where-    warned = unsafePerformIO $ newIORef False-    {-# NOINLINE warned #-} --- | Seed a PRNG with data from the system's fast source of pseudo-random--- numbers. All the caveats of 'withSystemRandom' apply here as well.+seedCreatetionWarned :: IORef Bool+seedCreatetionWarned = unsafePerformIO $ newIORef False+{-# NOINLINE seedCreatetionWarned #-}++++-- | Generate random seed for generator using system's fast source of+--   pseudo-random numbers.+--+-- @since 0.15.0.0+createSystemSeed :: IO Seed+createSystemSeed = do+  seed <- createSystemRandomList+  return $! toSeed $ I.fromList seed++-- | Seed a PRNG with data from the system's fast source of+--   pseudo-random numbers. createSystemRandom :: IO GenIO-createSystemRandom = withSystemRandom (return :: GenIO -> IO GenIO)+createSystemRandom = initialize . I.fromList =<< createSystemRandomList ++-- | Seed PRNG with data from the system's fast source of+--   pseudo-random numbers and execute computation in ST monad.+--+-- @since 0.15.0.0+withSystemRandomST :: (forall s. Gen s -> ST s a) -> IO a+withSystemRandomST act = do+  seed <- createSystemSeed+  return $! runST $ act =<< restore seed++-- | Seed a PRNG with data from the system's fast source of+--   pseudo-random numbers, then run the given action.+--+--   This function is unsafe and for example allows STRefs or any+--   other mutable data structure to escape scope:+--+--   >>> ref <- withSystemRandom $ \_ -> newSTRef 1+--   >>> withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref+--   2+--   >>> withSystemRandom $ \_ -> modifySTRef ref succ >> readSTRef ref+--   3+withSystemRandom :: PrimBase m+                 => (Gen (PrimState m) -> m a) -> IO a+withSystemRandom act = do+  seed <- createSystemSeed+  unsafePrimToIO $ act =<< restore seed+{-# DEPRECATED withSystemRandom "Use withSystemRandomST or createSystemSeed or createSystemRandom instead" #-}++ -- | Compute the next index into the state pool.  This is simply -- addition modulo 256. nextIndex :: Integral a => a -> Int@@ -453,37 +565,19 @@ {-# INLINE aa #-}  uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32+-- NOTE [Carry value] uniformWord32 (Gen q) = do   i  <- nextIndex `liftM` M.unsafeRead q ioff   c  <- fromIntegral `liftM` M.unsafeRead q coff   qi <- fromIntegral `liftM` M.unsafeRead q i   let t  = aa * qi + c-      -- The comments in this function are a proof that:-      --   "if the carry value is strictly smaller than the multiplicator,-      --    the next carry value is also strictly smaller than the multiplicator."-      -- Eventhough the proof is written in terms of the actual value of the multiplicator,-      --   it holds for any multiplicator value /not/ greater than maxBound 'Word32'-      ---      --    (In the code, the multiplicator is aa, the carry value is c,-      --     the next carry value is c''.)-      ---      -- So we'll assume that c < aa, and show that c'' < aa :-      ---      -- by definition, aa = 0x5BCF5AB2, qi <= 0xFFFFFFFF (because it is a 'Word32')-      -- hence aa*qi <= 0x5BCF5AB200000000 - 0x5BCF5AB2.-      ---      -- hence t  < 0x5BCF5AB200000000 (because t = aa * qi + c and c < 0x5BCF5AB2)-      -- hence t <= 0x5BCF5AB1FFFFFFFF       c' = fromIntegral (t `shiftR` 32)-      --       c' <         0x5BCF5AB1       x  = fromIntegral t + c'       (# x', c'' #)  | x < c'    = (# x + 1, c' + 1 #)                      | otherwise = (# x,     c' #)-      -- hence c'' <        0x5BCF5AB2,-      -- hence c'' < aa, which is what we wanted to prove.   M.unsafeWrite q i x'   M.unsafeWrite q ioff (fromIntegral i)-  M.unsafeWrite q coff (fromIntegral c'')+  M.unsafeWrite q coff c''   return x' {-# INLINE uniformWord32 #-} @@ -513,7 +607,7 @@   M.unsafeWrite q i x'   M.unsafeWrite q j y'   M.unsafeWrite q ioff (fromIntegral j)-  M.unsafeWrite q coff (fromIntegral d'')+  M.unsafeWrite q coff d''   return $! f x' y' {-# INLINE uniform2 #-} @@ -532,24 +626,8 @@ type instance Unsigned Word32 = Word32 type instance Unsigned Word64 = Word64 --- This is workaround for bug #25.------ GHC-7.6 has a bug (#8072) which results in calculation of wrong--- number of buckets in function `uniformRange'. Consequently uniformR--- generates values in wrong range.------ Bug only affects 32-bit systems and Int/Word data types. Word32--- works just fine. So we set Word32 as unsigned counterpart for Int--- and Word on 32-bit systems. It's done only for GHC-7.6 because--- other versions are unaffected by the bug and we expect that GHC may--- optimise code which uses Word better.-#if (WORD_SIZE_IN_BITS < 64) && (__GLASGOW_HASKELL__ == 706)-type instance Unsigned Int   = Word32-type instance Unsigned Word  = Word32-#else type instance Unsigned Int   = Word type instance Unsigned Word  = Word-#endif   -- Subtract two numbers under assumption that x>=y and store result in@@ -591,12 +669,25 @@ -- | Generate a vector of pseudo-random variates.  This is not -- necessarily faster than invoking 'uniform' repeatedly in a loop, -- but it may be more convenient to use in some situations.-uniformVector :: (PrimMonad m, Variate a, Vector v a)-             => Gen (PrimState m) -> Int -> m (v a)-uniformVector gen n = G.replicateM n (uniform gen)+uniformVector+  :: (PrimMonad m, Random.StatefulGen g m, Random.Uniform a, Vector v a)+  => g -> Int -> m (v a)+-- NOTE: We use in-place mutation in order to generate vector instead+--       of generateM because latter will go though intermediate list until+--       we're working in IO/ST monad+--+-- See: https://github.com/haskell/vector/issues/208 for details+uniformVector gen n = do+  mu <- GM.unsafeNew n+  let go !i | i < n     = Random.uniformM gen >>= GM.unsafeWrite mu i >> go (i+1)+            | otherwise = G.unsafeFreeze mu+  go 0 {-# INLINE uniformVector #-}  +-- This is default seed for the generator and used when no seed is+-- specified or seed is only partial. It's not known how it was+-- generated but it looks random enough defaultSeed :: I.Vector Word32 defaultSeed = I.fromList [   0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184,@@ -695,3 +786,40 @@ -- -- This is almost as compact as the original code that the compiler -- rejected.++++-- $setup+--+-- >>> import Control.Monad+-- >>> import Data.Word+-- >>> import Data.STRef+-- >>> :set -Wno-deprecations+++-- NOTE [Carry value]+-- ------------------+-- This is proof of statement:+--+-- > if the carry value is strictly smaller than the multiplicator,+-- > the next carry value is also strictly smaller than the multiplicator.+--+-- Eventhough the proof is written in terms of the actual value of the+-- multiplicator, it holds for any multiplicator value /not/ greater+-- than maxBound 'Word32'+--+--    (In the code, the multiplicator is aa, the carry value is c,+--     the next carry value is c''.)+--+-- So we'll assume that c < aa, and show that c'' < aa :+--+-- by definition, aa = 0x5BCF5AB2, qi <= 0xFFFFFFFF (because it is a 'Word32')+--+-- Then we get following:+--+--    aa*qi <= 0x5BCF5AB200000000 - 0x5BCF5AB2.+--    t     <  0x5BCF5AB200000000 (because t = aa * qi + c and c < 0x5BCF5AB2)+--    t     <= 0x5BCF5AB1FFFFFFFF+--    c'    <  0x5BCF5AB1+--    c''   <  0x5BCF5AB2,+--    c''   < aa, which is what we wanted to prove.

@@ -29,7 +29,6 @@   ) where  import Control.Arrow           (second,(***))-import Control.Monad.Primitive (PrimMonad(..))  import Data.Word import Data.Int@@ -41,13 +40,12 @@ import qualified Data.Vector                 as V import Data.Vector.Generic (Vector) import Numeric.SpecFunctions (logFactorial)+import System.Random.Stateful  import Prelude hiding ((++)) -import System.Random.MWC  - -- | A lookup table for arbitrary discrete distributions. It allows -- the generation of random variates in /O(1)/. Note that probability -- is quantized in units of @1/2^32@, and all distributions with@@ -77,11 +75,10 @@   -- | Generate a random value using a condensed table.-genFromTable :: (PrimMonad m, Vector v a) =>-                CondensedTable v a -> Gen (PrimState m) -> m a+genFromTable :: (StatefulGen g m, Vector v a) => CondensedTable v a -> g -> m a {-# INLINE genFromTable #-} genFromTable table gen = do-  w <- uniform gen+  w <- uniformM gen   return $ lookupTable table $ fromIntegral (w :: Word32)  lookupTable :: Vector v a => CondensedTable v a -> Word64 -> a

@@ -47,7 +47,7 @@ #endif import Data.Traversable (mapM) import Data.Word (Word32)-import System.Random.MWC (Gen, uniform, uniformR)+import System.Random.Stateful (StatefulGen(..),Uniform(..),UniformRange(..),uniformDoublePositive01M) import qualified Data.Vector.Unboxed         as I import qualified Data.Vector.Generic         as G import qualified Data.Vector.Generic.Mutable as M@@ -58,10 +58,10 @@  -- | Generate a normally distributed random variate with given mean -- and standard deviation.-normal :: PrimMonad m+normal :: StatefulGen g m        => Double                -- ^ Mean        -> Double                -- ^ Standard deviation-       -> Gen (PrimState m)+       -> g        -> m Double {-# INLINE normal #-} normal m s gen = do@@ -75,13 +75,13 @@ -- Compared to the ziggurat algorithm usually used, this is slower, -- but generates more independent variates that pass stringent tests -- of randomness.-standard :: PrimMonad m => Gen (PrimState m) -> m Double+standard :: StatefulGen g m => g -> m Double {-# INLINE standard #-} standard gen = loop   where     loop = do-      u  <- (subtract 1 . (*2)) `liftM` uniform gen-      ri <- uniform gen+      u  <- (subtract 1 . (*2)) `liftM` uniformDoublePositive01M gen+      ri <- uniformM gen       let i  = fromIntegral ((ri :: Word32) .&. 127)           bi = I.unsafeIndex blocks i           bj = I.unsafeIndex blocks (i+1)@@ -93,14 +93,14 @@                  xx = x * x                  d  = exp (-0.5 * (bi * bi - xx))                  e  = exp (-0.5 * (bj * bj - xx))-             c <- uniform gen+             c <- uniformDoublePositive01M gen              if e + c * (d - e) < 1                then return x                else loop     normalTail neg  = tailing       where tailing  = do-              x <- ((/rNorm) . log) `liftM` uniform gen-              y <- log              `liftM` uniform gen+              x <- ((/rNorm) . log) `liftM` uniformDoublePositive01M gen+              y <- log              `liftM` uniformDoublePositive01M gen               if y * (-2) < x * x                 then tailing                 else return $! if neg then x - rNorm else rNorm - x@@ -129,36 +129,36 @@   -- | Generate an exponentially distributed random variate.-exponential :: PrimMonad m+exponential :: StatefulGen g m             => Double            -- ^ Scale parameter-            -> Gen (PrimState m) -- ^ Generator+            -> g                 -- ^ Generator             -> m Double {-# INLINE exponential #-} exponential b gen = do-  x <- uniform gen+  x <- uniformDoublePositive01M gen   return $! - log x / b   -- | Generate truncated exponentially distributed random variate.-truncatedExp :: PrimMonad m+truncatedExp :: StatefulGen g m              => Double            -- ^ Scale parameter              -> (Double,Double)   -- ^ Range to which distribution is                                   --   truncated. Values may be negative.-             -> Gen (PrimState m) -- ^ Generator.+             -> g                 -- ^ Generator.              -> m Double {-# INLINE truncatedExp #-} truncatedExp scale (a,b) gen = do   -- We shift a to 0 and then generate distribution truncated to [0,b-a]   -- It's easier   let delta = b - a-  p <- uniform gen+  p <- uniformDoublePositive01M gen   return $! a - log ( (1 - p) + p*exp(-scale*delta)) / scale  -- | Random variate generator for gamma distribution.-gamma :: PrimMonad m+gamma :: (StatefulGen g m)       => Double                 -- ^ Shape parameter       -> Double                 -- ^ Scale parameter-      -> Gen (PrimState m)      -- ^ Generator+      -> g                      -- ^ Generator       -> m Double {-# INLINE gamma #-} gamma a b gen@@ -167,13 +167,13 @@     where       mainloop = do         T x v <- innerloop-        u     <- uniform gen+        u     <- uniformDoublePositive01M gen         let cont =  u > 1 - 0.331 * sqr (sqr x)                  && log u > 0.5 * sqr x + a1 * (1 - v + log v) -- Rarely evaluated         case () of           _| cont      -> mainloop            | a >= 1    -> return $! a1 * v * b-           | otherwise -> do y <- uniform gen+           | otherwise -> do y <- uniformDoublePositive01M gen                              return $! y ** (1 / a) * a1 * v * b       -- inner loop       innerloop = do@@ -188,9 +188,9 @@   -- | Random variate generator for the chi square distribution.-chiSquare :: PrimMonad m+chiSquare :: StatefulGen g m           => Int                -- ^ Number of degrees of freedom-          -> Gen (PrimState m)  -- ^ Generator+          -> g                  -- ^ Generator           -> m Double {-# INLINE chiSquare #-} chiSquare n gen@@ -201,14 +201,14 @@ -- | Random variate generator for the geometric distribution, -- computing the number of failures before success. Distribution's -- support is [0..].-geometric0 :: PrimMonad m+geometric0 :: StatefulGen g m            => Double            -- ^ /p/ success probability lies in (0,1]-           -> Gen (PrimState m) -- ^ Generator+           -> g                 -- ^ Generator            -> m Int {-# INLINE geometric0 #-} geometric0 p gen   | p == 1          = return 0-  | p >  0 && p < 1 = do q <- uniform gen+  | p >  0 && p < 1 = do q <- uniformDoublePositive01M gen                          -- FIXME: We want to use log1p here but it will                          --        introduce dependency on math-functions.                          return $! floor $ log q / log (1 - p)@@ -217,19 +217,19 @@ -- | Random variate generator for geometric distribution for number of -- trials. Distribution's support is [1..] (i.e. just 'geometric0' -- shifted by 1).-geometric1 :: PrimMonad m+geometric1 :: StatefulGen g m            => Double            -- ^ /p/ success probability lies in (0,1]-           -> Gen (PrimState m) -- ^ Generator+           -> g                 -- ^ Generator            -> m Int {-# INLINE geometric1 #-} geometric1 p gen = do n <- geometric0 p gen                       return $! n + 1  -- | Random variate generator for Beta distribution-beta :: PrimMonad m+beta :: StatefulGen g m      => Double            -- ^ alpha (>0)      -> Double            -- ^ beta  (>0)-     -> Gen (PrimState m) -- ^ Generator+     -> g                 -- ^ Generator      -> m Double {-# INLINE beta #-} beta a b gen = do@@ -238,9 +238,9 @@   return $! x / (x+y)  -- | Random variate generator for Dirichlet distribution-dirichlet :: (PrimMonad m, Traversable t)+dirichlet :: (StatefulGen g m, Traversable t)           => t Double          -- ^ container of parameters-          -> Gen (PrimState m) -- ^ Generator+          -> g                 -- ^ Generator           -> m (t Double) {-# INLINE dirichlet #-} dirichlet t gen = do@@ -249,12 +249,12 @@   return $ fmap (/total) t'  -- | Random variate generator for Bernoulli distribution-bernoulli :: PrimMonad m+bernoulli :: StatefulGen g m           => Double            -- ^ Probability of success (returning True)-          -> Gen (PrimState m) -- ^ Generator+          -> g                 -- ^ Generator           -> m Bool {-# INLINE bernoulli #-}-bernoulli p gen = (<p) `liftM` uniform gen+bernoulli p gen = (<p) `liftM` uniformDoublePositive01M gen  -- | Random variate generator for categorical distribution. --@@ -262,16 +262,16 @@ --   "System.Random.MWC.CondensedTable" will offer better --   performance.  If only few is needed this function will faster --   since it avoids costs of setting up table.-categorical :: (PrimMonad m, G.Vector v Double)+categorical :: (StatefulGen g m, G.Vector v Double)             => v Double          -- ^ List of weights [>0]-            -> Gen (PrimState m) -- ^ Generator+            -> g                 -- ^ Generator             -> m Int {-# INLINE categorical #-} categorical v gen     | G.null v = pkgError "categorical" "empty weights!"     | otherwise = do         let cv  = G.scanl1' (+) v-        p <- (G.last cv *) `liftM` uniform gen+        p <- (G.last cv *) `liftM` uniformDoublePositive01M gen         return $! case G.findIndex (>=p) cv of                     Just i  -> i                     Nothing -> pkgError "categorical" "bad weights!"@@ -279,9 +279,9 @@ -- | Random variate generator for categorical distribution where the --   weights are in the log domain. It's implemented in terms of --   'categorical'.-logCategorical :: (PrimMonad m, G.Vector v Double)+logCategorical :: (StatefulGen g m, G.Vector v Double)                => v Double          -- ^ List of logarithms of weights-               -> Gen (PrimState m) -- ^ Generator+               -> g                 -- ^ Generator                -> m Int {-# INLINE logCategorical #-} logCategorical v gen@@ -294,9 +294,9 @@ --   It returns random permutation of vector /[0 .. n-1]/. -- --   This is the Fisher-Yates shuffle-uniformPermutation :: forall m v. (PrimMonad m, G.Vector v Int)+uniformPermutation :: forall g m v. (StatefulGen g m, PrimMonad m, G.Vector v Int)                    => Int-                   -> Gen (PrimState m)+                   -> g                    -> m (v Int) {-# INLINE uniformPermutation #-} uniformPermutation n gen@@ -306,9 +306,9 @@ -- | Random variate generator for a uniformly distributed shuffle (all --   shuffles are equiprobable) of a vector. It uses Fisher-Yates --   shuffle algorithm.-uniformShuffle :: (PrimMonad m, G.Vector v a)+uniformShuffle :: (StatefulGen g m, PrimMonad m, G.Vector v a)                => v a-               -> Gen (PrimState m)+               -> g                -> m (v a) {-# INLINE uniformShuffle #-} uniformShuffle vec gen@@ -320,9 +320,9 @@  -- | In-place uniformly distributed shuffle (all shuffles are --   equiprobable)of a vector.-uniformShuffleM :: (PrimMonad m, M.MVector v a)+uniformShuffleM :: (StatefulGen g m, PrimMonad m, M.MVector v a)                 => v (PrimState m) a-                -> Gen (PrimState m)+                -> g                 -> m () {-# INLINE uniformShuffleM #-} uniformShuffleM vec gen@@ -332,7 +332,7 @@     n   = M.length vec     lst = n-1     loop i | i == lst  = return ()-           | otherwise = do j <- uniformR (i,lst) gen+           | otherwise = do j <- uniformRM (i,lst) gen                             M.unsafeSwap vec i j                             loop (i+1) @@ -344,8 +344,6 @@ pkgError :: String -> String -> a pkgError func msg = error $ "System.Random.MWC.Distributions." ++ func ++                             ": " ++ msg--  -- $references --

@@ -0,0 +1,111 @@+{-# LANGUAGE BangPatterns #-}+import Control.Exception+import Control.Monad+import Control.Monad.ST+import Gauge.Main+import Data.Int+import Data.Word+import qualified Data.Vector.Unboxed as U+import qualified System.Random as R+import System.Random.MWC+import System.Random.MWC.Distributions+import System.Random.MWC.CondensedTable+import qualified System.Random.Mersenne as M++makeTableUniform :: Int -> CondensedTable U.Vector Int+makeTableUniform n =+  tableFromProbabilities $ U.zip (U.enumFromN 0 n) (U.replicate n (1 / fromIntegral n))+{-# INLINE makeTableUniform #-}+++main = do+  mwc <- create+  mtg <- M.newMTGen . Just =<< uniform mwc+  defaultMain+    [ bgroup "mwc"+      -- One letter group names are used so they will fit on the plot.+      --+      --  U - uniform+      --  R - uniformR+      --  D - distribution+      [ bgroup "U"+        [ bench "Double"  $ nfIO (uniform mwc :: IO Double)+        , bench "Int"     $ nfIO (uniform mwc :: IO Int)+        , bench "Int8"    $ nfIO (uniform mwc :: IO Int8)+        , bench "Int16"   $ nfIO (uniform mwc :: IO Int16)+        , bench "Int32"   $ nfIO (uniform mwc :: IO Int32)+        , bench "Int64"   $ nfIO (uniform mwc :: IO Int64)+        , bench "Word"    $ nfIO (uniform mwc :: IO Word)+        , bench "Word8"   $ nfIO (uniform mwc :: IO Word8)+        , bench "Word16"  $ nfIO (uniform mwc :: IO Word16)+        , bench "Word32"  $ nfIO (uniform mwc :: IO Word32)+        , bench "Word64"  $ nfIO (uniform mwc :: IO Word64)+        ]+      , bgroup "R"+        -- I'm not entirely convinced that this is right way to test+        -- uniformR. /A.Khudyakov/+        [ bench "Double"  $ nfIO (uniformR (-3.21,26) mwc :: IO Double)+        , bench "Int"     $ nfIO (uniformR (-12,679)  mwc :: IO Int)+        , bench "Int8"    $ nfIO (uniformR (-12,4)    mwc :: IO Int8)+        , bench "Int16"   $ nfIO (uniformR (-12,679)  mwc :: IO Int16)+        , bench "Int32"   $ nfIO (uniformR (-12,679)  mwc :: IO Int32)+        , bench "Int64"   $ nfIO (uniformR (-12,679)  mwc :: IO Int64)+        , bench "Word"    $ nfIO (uniformR (34,633)   mwc :: IO Word)+        , bench "Word8"   $ nfIO (uniformR (34,63)    mwc :: IO Word8)+        , bench "Word16"  $ nfIO (uniformR (34,633)   mwc :: IO Word16)+        , bench "Word32"  $ nfIO (uniformR (34,633)   mwc :: IO Word32)+        , bench "Word64"  $ nfIO (uniformR (34,633)   mwc :: IO Word64)+        ]+      , bgroup "D"+        [ bench "standard"    $ nfIO (standard      mwc :: IO Double)+        , bench "normal"      $ nfIO (normal 1 3    mwc :: IO Double)+          -- Regression tests for #16. These functions should take 10x+          -- longer to execute.+          --+          -- N.B. Bang patterns are necessary to trigger the bug with+          --      GHC 7.6+        , bench "standard/N"  $ nfIO $ replicateM_ 10 $ do+                                 !_ <- standard mwc :: IO Double+                                 return ()+        , bench "normal/N"    $ nfIO $ replicateM_ 10 $ do+                                 !_ <- normal 1 3 mwc :: IO Double+                                 return ()+        , bench "exponential" $ nfIO (exponential 3 mwc :: IO Double)+        , bench "gamma,a<1"   $ nfIO (gamma 0.5 1   mwc :: IO Double)+        , bench "gamma,a>1"   $ nfIO (gamma 2   1   mwc :: IO Double)+        , bench "chiSquare"   $ nfIO (chiSquare 4   mwc :: IO Double)+        ]+      , bgroup "CT/gen" $ concat+        [ [ bench ("uniform "++show i) $ nfIO (genFromTable (makeTableUniform i) mwc :: IO Int)+          | i <- [2..10]+          ]+        , [ bench ("poisson " ++ show l) $ nfIO (genFromTable (tablePoisson l) mwc :: IO Int)+          | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]+          ]+        , [ bench ("binomial " ++ show p ++ " " ++ show n) $ nfIO (genFromTable (tableBinomial n p) mwc :: IO Int)+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]+          ]+        ]+      , bgroup "CT/table" $ concat+        [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i+          | i <- [2..30]+          ]+        , [ bench ("poisson " ++ show l) $ whnf tablePoisson l+          | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]+          ]+        , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnf (tableBinomial n) p+          | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]+          ]+        ]+      ]+    , bgroup "random"+      [+        bench "Double" $ nfIO (R.randomIO >>= evaluate :: IO Double)+      , bench "Int"    $ nfIO (R.randomIO >>= evaluate :: IO Int)+      ]+    , bgroup "mersenne"+      [+        bench "Double" $ nfIO (M.random mtg :: IO Double)+      , bench "Int"    $ nfIO (M.random mtg :: IO Int)+      ]+    ]

@@ -1,3 +1,22 @@+## Changes in 0.15.0.0++  * `withSystemRandomST` and `createSystemSeed` are added.++  * `withSystemRandom` is deprecated.++  * `random>=1.2` is dependency of `mwc-random`.++  * Instances for type classes `StatefulGen` & `FrozenGen` defined in random-1.2+    are added for `Gen`.++  * Functions in `System.Random.MWC.Distributions` and+    `System.Random.MWC.CondensedTable` now work with arbitrary `StatefulGen`++  * `System.Random.MWC.uniformVector` now works with arbitrary `StatefulGen` as+    well and uses in-place initialization instead of `generateM`. It should be+    faster for anything but IO and ST (those shoud remain same).++ ## Changes in 0.14.0.0    * Low level functions for acquiring random data for initialization

@@ -1,5 +1,5 @@ name:           mwc-random-version:        0.14.0.0+version:        0.15.0.0 synopsis:       Fast, high quality pseudo random number generation description:   This package contains code for generating high quality random@@ -24,18 +24,30 @@ copyright:      2009, 2010, 2011 Bryan O'Sullivan category:       Math, Statistics build-type:     Simple-cabal-version:  >= 1.8.0.4+cabal-version:  >= 1.10 extra-source-files:   changelog.md   README.markdown +tested-with:+    GHC ==7.10.3+     || ==8.0.2+     || ==8.2.2+     || ==8.4.4+     || ==8.6.5+     || ==8.8.3+     || ==8.10.1+  , GHCJS ==8.4+ library+  default-language: Haskell2010   exposed-modules: System.Random.MWC                    System.Random.MWC.Distributions                    System.Random.MWC.CondensedTable                    System.Random.MWC.SeedSource-  build-depends: base           >= 4.5 && < 5+  build-depends: base           >= 4.8 && < 5                , primitive      >= 0.6+               , random         >= 1.2                , time                , vector         >= 0.7                , math-functions >= 0.2.1.0@@ -50,3 +62,50 @@ source-repository head   type:     mercurial   location: https://bitbucket.org/bos/mwc-random+++benchmark mwc-bench+  type:           exitcode-stdio-1.0+  hs-source-dirs: bench+  main-is:        Benchmark.hs+  default-language: Haskell2010+  build-depends: base < 5+               , vector          >= 0.11+               , gauge           >= 0.2.5+               , mersenne-random+               , mwc-random+               , random+++test-suite mwc-prop-tests+  type:           exitcode-stdio-1.0+  hs-source-dirs: tests+  main-is:        props.hs+  default-language: Haskell2010+  ghc-options:+    -Wall -threaded -rtsopts++  build-depends: base+               , mwc-random+               , QuickCheck                 >=2.2+               , vector                     >=0.12.1+               , tasty >=1.3.1+               , tasty-quickcheck+               , tasty-hunit++test-suite mwc-doctests+  type:             exitcode-stdio-1.0+  main-is:          doctests.hs+  hs-source-dirs:   tests+  default-language: Haskell2010+  if impl(ghcjs) || impl(ghc < 8.0)+    Buildable: False+  build-depends:+            base       -any+          , mwc-random -any+          , doctest    >=0.15 && <0.18+            --+          , bytestring+          , primitive+          , vector     >=0.11+          , random     >=1.2

@@ -0,0 +1,4 @@+import Test.DocTest (doctest)++main :: IO ()+main = doctest ["-fobject-code", "System/Random/MWC.hs"]

@@ -0,0 +1,120 @@+import Control.Monad+import Data.Word+import qualified Data.Vector.Unboxed as U++import Test.Tasty+import Test.Tasty.QuickCheck+import Test.Tasty.HUnit+import Test.QuickCheck.Monadic++import System.Random.MWC+++----------------------------------------------------------------+--+----------------------------------------------------------------++main :: IO ()+main = do+  g0 <- createSystemRandom+  defaultMain $ testGroup "mwc"+    [ testProperty "save/restore"      $ prop_SeedSaveRestore g0+    , testCase     "user save/restore" $ saveRestoreUserSeed+    , testCase     "empty seed data"   $ emptySeed+    , testCase     "output correct"    $ do+        g  <- create+        xs <- replicateM 513 (uniform g)+        assertEqual "[Word32]" xs golden+    ]++updateGenState :: GenIO -> IO ()+updateGenState g = replicateM_ 256 (uniform g :: IO Word32)+++prop_SeedSaveRestore :: GenIO -> Property+prop_SeedSaveRestore g = monadicIO  $ do+  run $ updateGenState g+  seed  <- run $ save g+  seed' <- run $ save =<< restore seed+  return $ seed == seed'++saveRestoreUserSeed :: IO ()+saveRestoreUserSeed = do+  let seed = toSeed $ U.replicate 258 0+  seed' <- save =<< restore seed+  assertEqual "Seeds must be equal" seed' seed+      +emptySeed :: IO ()+emptySeed = do+  let seed = toSeed U.empty+  seed' <- save =<< create+  assertEqual "Seeds must be equal" seed' seed++-- First 513 values generated from seed made using create+golden :: [Word32]+golden =+  [ 2254043345, 562229898, 1034503294, 2470032534, 2831944869, 3042560015, 838672965, 715056843+  , 3122641307, 2300516242, 4079538318, 3722020688, 98524204, 1450170923, 2669500465, 2890402829+  , 114212910, 1914313000, 2389251496, 116282477, 1771812561, 1606473512, 1086420617, 3652430775+  , 1165083752, 3599954795, 3006722175, 341614641, 3000394300, 1378097585, 1551512487, 81211762+  , 604209599, 3949866361, 77745071, 3170410267, 752447516, 1213023833, 1624321744, 3251868348+  , 1584957570, 2296897736, 3305840056, 1158966242, 2458014362, 1919777052, 3203159823, 3230279656+  , 755741068, 3005087942, 2478156967, 410224731, 1196248614, 3302310440, 3295868805, 108051054+  , 1010042411, 2725695484, 2201528637, 667561409, 79601486, 50029770, 566202616, 3217300833+  , 2162817014, 925506837, 1527015413, 3079491438, 927252446, 118306579, 499811870, 2973454232+  , 2979271640, 4078978924, 1864075883, 197741457, 296365782, 1784247291, 236572186, 464208268+  , 1769568958, 827682258, 4247376295, 2959098022, 1183860331, 2475064236, 3952901213, 1953014945+  , 393081236, 1616500498, 2201176136, 1663813362, 2167124739, 630903810, 113470040, 924745892+  , 1081531735, 4039388931, 4118728223, 107819176, 2212875141, 1941653033, 3660517172, 192973521+  , 3653156164, 1878601439, 3028195526, 2545631291, 3882334975, 456082861, 2775938704, 3813508885+  , 1758481462, 3332769695, 3595846251, 3745876981, 152488869, 2555728588, 3058747945, 39382408+  , 520595021, 2185388418, 3502636573, 2650173199, 1077668433, 3548643646, 71562049, 2726649517+  , 494210825, 1208915815, 620990806, 2877290965, 3253243521, 804166732, 2481889113, 623399529+  , 44880343, 183645859, 3283683418, 2214754452, 419328482, 4224066437, 1102669380, 1997964721+  , 2437245376, 985749802, 858381069, 116806511, 1771295365, 97352549, 341972923, 2971905841+  , 110707773, 950500868, 1237119233, 691764764, 896381812, 1528998276, 1269357470, 2567094423+  , 52141189, 2722993417, 80628658, 3919817965, 3615946076, 899371181, 46940285, 4010779728+  , 318101834, 30736609, 3577200709, 971882724, 1478800972, 3769640027, 3706909300, 3300631811+  , 4057825972, 4285058790, 2329759553, 2967563409, 4080096760, 2762613004, 2518395275, 295718526+  , 598435593, 2385852565, 2608425408, 604857293, 2246982455, 919156819, 1721573814, 2502545603+  , 643962859, 587823425, 3508582012, 1777595823, 4119929334, 2833342174, 414044876, 2469473258+  , 289159600, 3715175415, 966867024, 788102818, 3197534326, 3571396978, 3508903890, 570753009+  , 4273926277, 3301521986, 1411959102, 2766249515, 4071012597, 959442028, 1962463990, 1098904190+  , 714719899, 562204808, 1658783410, 1471669042, 2565780129, 1616648894, 4236521717, 1788863789+  , 3068674883, 191936470, 253084644, 1915647866, 276372665, 2117183118, 3704675319, 218791054+  , 3680045802, 406662689, 3844864229, 91140313, 3834015630, 25116147, 904830493, 3152559113+  , 820358622, 1301896358, 296152699, 2202014455, 4256659428, 1175171414, 3287520873, 2028006499+  , 327448717, 2095642873, 3798661296, 58567008, 3907537112, 3691259011, 1730142328, 2373011713+  , 3387040741, 3189417655, 2949233059, 1238379614, 1813238023, 1064726446, 1339055235, 1744523609+  , 279811576, 2934103599, 283542302, 994488448, 418691747, 1062780152, 102211875, 4071713296+  , 1790834038, 1035092527, 2374272359, 3558280982, 1927663822, 3645417844, 3481790745, 3566282546+  , 2000290859, 505518126, 363501589, 4075468679, 3247300709, 3705242654, 2731103609, 2836871038+  , 589640144, 2546495106, 84767518, 1376911639, 2400770705, 527489676, 3804134352, 150084021+  , 240070593, 3807594859, 3518576690, 659503830, 2239678479, 1273668921, 4271050554, 3090482972+  , 401956859, 1772128561, 4438455, 1989666158, 2521484677, 3960178700, 4220196277, 1033999035+  , 2214785840, 3428469341, 428564336, 2517446784, 3935757188, 3294001677, 1037971963, 3590324170+  , 1220969729, 1719719817, 807688972, 77076422, 4251553858, 3963852375, 326128795, 3277818295+  , 3671513069, 549617771, 1683950556, 3352913781, 409318429, 2456264774, 4036950639, 1162718475+  , 83888874, 5578966, 172866494, 1542278848, 455546979, 1296511553, 4263636440, 2450589064+  , 372411483, 211216338, 2632256495, 2393754408, 1336054289, 4087203071, 3159642437, 1933346856+  , 2914152714, 3805541979, 2769740793, 1161287028, 2289749561, 4124509890, 2128452935, 210531695+  , 4250709834, 390950534, 1421430300, 3030519715, 3228987297, 3086837053, 2866915453, 2335948692+  , 1684378991, 2575634059, 4153427304, 2426048796, 4197556954, 2605152326, 2909410733, 2424889219+  , 654577921, 811955499, 118126602, 504071559, 1278756230, 3896458168, 4105558075, 750276169+  , 1120805572, 1762689330, 993728154, 1104363215, 774344996, 4077568952, 2183487324, 994724370+  , 3323036885, 3880704963, 746305447, 961608310, 2030117337, 453935768, 800490463, 1034636+  , 2323633564, 602565693, 806061242, 1899269713, 162686347, 467541008, 1529175313, 282891502+  , 2529616339, 2930657178, 464272784, 2878535316, 807165854, 3209080518, 4080120278, 347748171+  , 3972126063, 284174728, 2498328933, 1723872460, 143845955, 4223866687, 1761495357, 1544646770+  , 4206103283, 3771574626, 642165282, 1119501013, 3514063332, 1443320304, 4056369796, 3602131475+  , 1422908288, 804093687, 431176780, 40108717, 2998264213, 3705835674, 169805085, 454593842+  , 2781536994, 2385225212, 4137367775, 2631435125, 2347082354, 629238010, 3283635219, 3815791831+  , 1340400558, 4061846985, 3803921868, 3196119096, 718610843, 3694290834, 2169960411, 2407155570+  , 2557480499, 16164105, 480957288, 2155919829, 2490067282, 2356287132, 511737296, 1602800634+  , 1802275249, 3316832299, 50286484, 2106622541, 2352302834, 2538374315, 344766394, 2777260569+  , 1215135803, 2229011963, 114632277, 1645499402, 1111617833, 3833259754, 928611385, 686744723+  , 1898396834, 2461932251, 2665457318, 3797019621, 868313114, 2366635205, 481934875, 1170532970+  , 642610859, 3150733309, 3508548582, 666714469, 711663449, 2436617656, 2681476315, 1637296693+  , 2487349478, 4174144946, 2793869557, 559398604, 1898140528, 991962870, 864792875, 3861665129+  , 4024051364, 3383200293, 773730975, 33517291, 2660126073, 689133464, 2248134097, 3874737781+  , 3358012678]