## 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))
(loop (cdr rs) (+ up 1) down))
(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.

Pages: 1 2

### 4 Responses to “RNG147”

1. Rutger said

In python 3

```def rng147(seed):
while True:
yield seed
seed = (147.0 * seed) % 1.0

r = rng147(0.1234)

for v in range(10):
print(next(r))
```
2. Graham said

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:

```module Main where

import           Data.Fixed  (mod')
import           Data.Stream (Stream)
import qualified Data.Stream as S

rng147 :: Double -> Stream Double
rng147 = S.iterate (\x -> (147 * x) `mod'` 1)

main :: IO ()
main = mapM_ print . S.take 10 . rng147 \$ pi / 10
```
3. Globules said

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 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))
```
```\$ ./rng147 10
0.14813489999999874
0.7758302999998143
4.705409997271204e-2
0.9169526959886696
0.7920463103344275
0.430807619160845
0.3287200166442119
0.32184244669915074
0.310839664775159
0.6934307219483742
```

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
\$
```
4. […] 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 […]