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.

Advertisements

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 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))
    
    $ ./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 […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: