RNG147
March 28, 2017
Let’s begin with a simple implementation of the random number generator; we’ll change this in a few moments:
(define rand
(let ((seed 0.3141592654))
(lambda args
(set! seed
(if (pair? args)
(car args)
(let ((x (* seed 147.0)))
(- x (floor x)))))
seed)))
Our random number generator returns the next random number in the sequence; if it is given an argument, it resets the seed:
> (rand) 0.18141201379999927 > (rand) 0.667566028599893 > (rand) 0.13220620418427131 > (rand 0.5) 0.5 > (rand) 0.5 > (rand) 0.5
Our random number works fine initially, but setting the seed to 0.5 causes it to loop repeatedly, as 147.0 * 0.5 = 73.5 and the fractional part of the calculation never changes. We can fix that by ensuring that simple seeds become complex:
(define rand
(let ((seed 0.3141592654))
(lambda args
(set! seed
(if (pair? args)
(sin (car args))
(let ((x (* seed 147.0)))
(- x (floor x)))))
seed)))
Taking the sine of the seed breaks the pattern:
> (rand) 0.18141201379999927 > (rand) 0.667566028599893 > (rand) 0.13220620418427131 > (rand 0.5) 0.479425538604203 > (rand) 0.47555417481784445 > (rand) 0.9064636982231349
That looks better. We’ll use that version of the random number generator for the remainder of the exercise. It is convenient to have a function that returns a list of random numbers; randlist returns a list of n random numbers, after optionally setting the seed of the first random number in the list (which is returned in reverse order, not that it matters):
(define (randlist n . args)
(when (pair? args) (rand (car args)))
(do ((n n (- n 1))
(rs (list) (cons (rand) rs)))
((zero? n) rs))))
We perform two tests to determine the suitability of the random number generator. The first test counts the frequency of the first k digits of the fraction:
(define (freq n k)
(uniq-c =
(sort < (map inexact->exact
(map floor
(map (lambda (x) (* x (expt 10 k)))
(randlist n)))))))
Here are results for k = 1 and k = 2:
> (freq 10000 1) ((0 . 983) (1 . 1067) (2 . 1019) (3 . 1040) (4 . 1043) (5 . 965) (6 . 994) (7 . 969) (8 . 928) (9 . 992)) > (freq 100000 2) ((0 . 981) (1 . 1040) (2 . 965) (3 . 1011) (4 . 930) (5 . 967) (6 . 1007) (7 . 987) (8 . 962) (9 . 966) (10 . 1036) (11 . 1028) (12 . 1026) (13 . 986) (14 . 958) (15 . 960) (16 . 999) (17 . 1048) (18 . 1067) (19 . 968) (20 . 1046) (21 . 930) (22 . 1014) (23 . 980) (24 . 952) (25 . 968) (26 . 1036) (27 . 960) (28 . 969) (29 . 1002) (30 . 998) (31 . 1006) (32 . 1001) (33 . 1074) (34 . 958) (35 . 1016) (36 . 988) (37 . 968) (38 . 994) (39 . 993) (40 . 981) (41 . 1094) (42 . 985) (43 . 999) (44 . 1010) (45 . 1003) (46 . 1016) (47 . 1017) (48 . 1008) (49 . 988) (50 . 942) (51 . 971) (52 . 1030) (53 . 1034) (54 . 1042) (55 . 997) (56 . 959) (57 . 1038) (58 . 1033) (59 . 1009) (60 . 1001) (61 . 1066) (62 . 998) (63 . 920) (64 . 983) (65 . 952) (66 . 982) (67 . 1014) (68 . 1003) (69 . 1021) (70 . 933) (71 . 1057) (72 . 1043) (73 . 1063) (74 . 1034) (75 . 1044) (76 . 966) (77 . 1039) (78 . 1034) (79 . 1020) (80 . 984) (81 . 970) (82 . 981) (83 . 1022) (84 . 1057) (85 . 983) (86 . 1024) (87 . 951) (88 . 973) (89 . 1009) (90 . 1017) (91 . 970) (92 . 991) (93 . 1020) (94 . 1003) (95 . 945) (96 . 988) (97 . 996) (98 . 1021) (99 . 1021))
We expect all of the counts to be about 1000, and they are. The second test looks at whether the sequence increases or decreases at each step:
(define (dir n)
(let loop ((rs (randlist n)) (up 0) (down 0))
(cond ((null? (cdr rs))
(list up (- n up down 1) down))
((< (car rs) (cadr rs))
(loop (cdr rs) (+ up 1) down))
((< (cadr rs) (car rs))
(loop (cdr rs) up (+ down 1)))
(else (loop (cdr rs) up down)))))
The test shows increases and decreases are about the same, and there are no ties:
> (dir 100000) (50036 0 49963)
There are many other tests that random number generators are typically required to pass, but, at first glance, RNG147 looks suitable for basic use in simulations that require random floating-point numbers, and is simple to implement and fast to operate. You can run the program at http://ideone.com/AUVQPw.
Programming Praxis 
In python 3
Your sine trick only moves the problem elsewhere; try seeding with pi/6 (arcsin of 1/2).
Solution in Haskell, using the Stream package for explicitly infinite lists:
Here’s another Haskell version. With an integer argument, n, it prints the first n random numbers. With no arguments it converts the stream of random Doubles into a stream of 32-bit unsigned integers. The purpose is to feed that binary data to the dieharder program, which evaluates its randomness. We don’t make use of all the bits in the Double; it’s just an easy way of getting the data to dieharder in a form that it likes.
{-# LANGUAGE ScopedTypeVariables #-} import Control.Monad (void) import Data.Word (Word32) import Foreign.Marshal.Alloc (alloca) import Foreign.Storable (poke, sizeOf) import System.Environment (getArgs) import System.IO (hPutBuf, stdout) -- A stream of "random" numbers. (We take the tail to exclude the seed from the -- values.) rng147 :: RealFrac b => b -> [b] rng147 = tail . iterate step where step x = let (_ :: Int, f) = properFraction (147 * x) in f main :: IO () main = do let rs = rng147 (0.1234567 :: Double) ns <- fmap (map read) getArgs :: IO [Int] case ns of -- Output an endless stream of 32-bit binary data for dieharder. [] -> put $ map cvt rs -- Print the first n random values. [n] -> mapM_ print $ take n rs -- Invalid argument(s). _ -> putStrLn "Usage: rng147 [n]" -- Convert a value to an 32-bit unsigned integer, possibly throwing away the -- least significant bits of the argument. cvt :: RealFrac a => a -> Word32 cvt x = round $ fromIntegral (maxBound :: Word32) * x -- Output the unsigned 32-bit values as binary data. put :: [Word32] -> IO () put is = let n = sizeOf (0 :: Word32) in alloca (\p -> void (mapM_ (\i -> poke p i >> hPutBuf stdout p n) is))Here’s the result of feeding the output to dieharder. I won’t try to interpret the results, here. (The -a argument tells it to run all the tests it has; -g 200 says that stdin is a series of 32-bit unsigned integers.) The summary is that 82 tests passed, 27 failed and 5 were weak.
$ ./rng147 | dieharder -a -g 200 #=============================================================================# # dieharder version 3.31.1 Copyright 2003 Robert G. Brown # #=============================================================================# rng_name |rands/second| Seed | stdin_input_raw| 1.01e+07 |2776686247| #=============================================================================# test_name |ntup| tsamples |psamples| p-value |Assessment #=============================================================================# diehard_birthdays| 0| 100| 100|0.64044752| PASSED diehard_operm5| 0| 1000000| 100|0.00000000| FAILED diehard_rank_32x32| 0| 40000| 100|0.45366134| PASSED diehard_rank_6x8| 0| 100000| 100|0.97658733| PASSED diehard_bitstream| 0| 2097152| 100|0.72508491| PASSED diehard_opso| 0| 2097152| 100|0.00000000| FAILED diehard_oqso| 0| 2097152| 100|0.00000000| FAILED diehard_dna| 0| 2097152| 100|0.00000000| FAILED diehard_count_1s_str| 0| 256000| 100|0.00000000| FAILED diehard_count_1s_byt| 0| 256000| 100|0.00000000| FAILED diehard_parking_lot| 0| 12000| 100|0.00000000| FAILED diehard_2dsphere| 2| 8000| 100|0.00000000| FAILED diehard_3dsphere| 3| 4000| 100|0.00000000| FAILED diehard_squeeze| 0| 100000| 100|0.00000000| FAILED diehard_sums| 0| 100| 100|0.06709806| PASSED diehard_runs| 0| 100000| 100|0.00000000| FAILED diehard_runs| 0| 100000| 100|0.00000000| FAILED diehard_craps| 0| 200000| 100|0.13428702| PASSED diehard_craps| 0| 200000| 100|0.00096831| WEAK marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED marsaglia_tsang_gcd| 0| 10000000| 100|0.00000000| FAILED sts_monobit| 1| 100000| 100|0.03994929| PASSED sts_runs| 2| 100000| 100|0.22979976| PASSED sts_serial| 1| 100000| 100|0.98471576| PASSED sts_serial| 2| 100000| 100|0.31624292| PASSED sts_serial| 3| 100000| 100|0.13274265| PASSED sts_serial| 3| 100000| 100|0.06274970| PASSED sts_serial| 4| 100000| 100|0.10906298| PASSED sts_serial| 4| 100000| 100|0.71359846| PASSED sts_serial| 5| 100000| 100|0.00790811| PASSED sts_serial| 5| 100000| 100|0.07231695| PASSED sts_serial| 6| 100000| 100|0.17659272| PASSED sts_serial| 6| 100000| 100|0.18513720| PASSED sts_serial| 7| 100000| 100|0.04315090| PASSED sts_serial| 7| 100000| 100|0.45634930| PASSED sts_serial| 8| 100000| 100|0.14826946| PASSED sts_serial| 8| 100000| 100|0.42520026| PASSED sts_serial| 9| 100000| 100|0.71305937| PASSED sts_serial| 9| 100000| 100|0.20822407| PASSED sts_serial| 10| 100000| 100|0.97154428| PASSED sts_serial| 10| 100000| 100|0.69812982| PASSED sts_serial| 11| 100000| 100|0.53454147| PASSED sts_serial| 11| 100000| 100|0.08444014| PASSED sts_serial| 12| 100000| 100|0.22453592| PASSED sts_serial| 12| 100000| 100|0.72236320| PASSED sts_serial| 13| 100000| 100|0.78366045| PASSED sts_serial| 13| 100000| 100|0.37616667| PASSED sts_serial| 14| 100000| 100|0.95764809| PASSED sts_serial| 14| 100000| 100|0.14472891| PASSED sts_serial| 15| 100000| 100|0.43901688| PASSED sts_serial| 15| 100000| 100|0.76714556| PASSED sts_serial| 16| 100000| 100|0.58587727| PASSED sts_serial| 16| 100000| 100|0.67425978| PASSED rgb_bitdist| 1| 100000| 100|0.05846675| PASSED rgb_bitdist| 2| 100000| 100|0.00004623| WEAK rgb_bitdist| 3| 100000| 100|0.66589418| PASSED rgb_bitdist| 4| 100000| 100|0.00026565| WEAK rgb_bitdist| 5| 100000| 100|0.07969442| PASSED rgb_bitdist| 6| 100000| 100|0.17850326| PASSED rgb_bitdist| 7| 100000| 100|0.25911992| PASSED rgb_bitdist| 8| 100000| 100|0.00045592| WEAK rgb_bitdist| 9| 100000| 100|0.00642740| PASSED rgb_bitdist| 10| 100000| 100|0.27514527| PASSED rgb_bitdist| 11| 100000| 100|0.91420132| PASSED rgb_bitdist| 12| 100000| 100|0.09738292| PASSED rgb_minimum_distance| 2| 10000| 1000|0.00000000| FAILED rgb_minimum_distance| 3| 10000| 1000|0.00000000| FAILED rgb_minimum_distance| 4| 10000| 1000|0.00000000| FAILED rgb_minimum_distance| 5| 10000| 1000|0.01232138| PASSED rgb_permutations| 2| 100000| 100|0.55323173| PASSED rgb_permutations| 3| 100000| 100|0.00000000| FAILED rgb_permutations| 4| 100000| 100|0.00000000| FAILED rgb_permutations| 5| 100000| 100|0.00000000| FAILED rgb_lagged_sum| 0| 1000000| 100|0.76792249| PASSED rgb_lagged_sum| 1| 1000000| 100|0.59838659| PASSED rgb_lagged_sum| 2| 1000000| 100|0.95607401| PASSED rgb_lagged_sum| 3| 1000000| 100|0.28565715| PASSED rgb_lagged_sum| 4| 1000000| 100|0.95085961| PASSED rgb_lagged_sum| 5| 1000000| 100|0.45582064| PASSED rgb_lagged_sum| 6| 1000000| 100|0.60565985| PASSED rgb_lagged_sum| 7| 1000000| 100|0.04195406| PASSED rgb_lagged_sum| 8| 1000000| 100|0.30831976| PASSED rgb_lagged_sum| 9| 1000000| 100|0.36546234| PASSED rgb_lagged_sum| 10| 1000000| 100|0.19936994| PASSED rgb_lagged_sum| 11| 1000000| 100|0.07303233| PASSED rgb_lagged_sum| 12| 1000000| 100|0.34645127| PASSED rgb_lagged_sum| 13| 1000000| 100|0.29195838| PASSED rgb_lagged_sum| 14| 1000000| 100|0.30094975| PASSED rgb_lagged_sum| 15| 1000000| 100|0.99790005| WEAK rgb_lagged_sum| 16| 1000000| 100|0.75196523| PASSED rgb_lagged_sum| 17| 1000000| 100|0.38809626| PASSED rgb_lagged_sum| 18| 1000000| 100|0.61180091| PASSED rgb_lagged_sum| 19| 1000000| 100|0.19440479| PASSED rgb_lagged_sum| 20| 1000000| 100|0.57917094| PASSED rgb_lagged_sum| 21| 1000000| 100|0.83129601| PASSED rgb_lagged_sum| 22| 1000000| 100|0.57292407| PASSED rgb_lagged_sum| 23| 1000000| 100|0.10381929| PASSED rgb_lagged_sum| 24| 1000000| 100|0.93651195| PASSED rgb_lagged_sum| 25| 1000000| 100|0.50848729| PASSED rgb_lagged_sum| 26| 1000000| 100|0.81756629| PASSED rgb_lagged_sum| 27| 1000000| 100|0.29194323| PASSED rgb_lagged_sum| 28| 1000000| 100|0.14934239| PASSED rgb_lagged_sum| 29| 1000000| 100|0.96722837| PASSED rgb_lagged_sum| 30| 1000000| 100|0.41045537| PASSED rgb_lagged_sum| 31| 1000000| 100|0.33379469| PASSED rgb_lagged_sum| 32| 1000000| 100|0.67528985| PASSED rgb_kstest_test| 0| 10000| 1000|0.12422629| PASSED dab_bytedistrib| 0| 51200000| 1|1.00000000| FAILED dab_dct| 256| 50000| 1|0.00000000| FAILED Preparing to run test 207. ntuple = 0 dab_filltree| 32| 15000000| 1|0.00000000| FAILED dab_filltree| 32| 15000000| 1|0.00000000| FAILED Preparing to run test 208. ntuple = 0 dab_filltree2| 0| 5000000| 1|0.00000000| FAILED dab_filltree2| 1| 5000000| 1|0.00000000| FAILED Preparing to run test 209. ntuple = 0 dab_monobit2| 12| 65000000| 1|1.00000000| FAILED $[…] built several random number generators: [1], [2], [3], [4], [5], [6], [7], [8], [9] (I didn’t realize it was so many until I went back and looked). In today’s exercise we […]