Back in the 1940s the mathematician John von Neumann invented the ‘middle square’ method of generating pseudo random numbers. The algorithm starts with a 4 digit number which is then squared. If the resulting number has fewer than 8 digits then it is padded with leading zeros. From this 8 digit number the middle 4 digits are extracted and returned as the result and also as input to repeating the algorithm.
This implementation, in Haskell, will manipulate the digits of the number as a list and so in outline we need to take the seed number and
‘square it‘ then ‘convert number to digits‘ then ‘pad to 8‘ then ‘get middle 4‘ then ‘convert digits to number‘
i.e. we can write a von Neumann Random number generator vnRan as
|
1 2 3 |
-- vnRan :: Int -> Int vnRan = digsToNum . mid4 . pad . numTodigs . sq |
Starting from the right… sq can be defined as (^2)
|
1 2 3 |
-- sq :: Int -> Int sq = (^ 2) |
numToDigs combines div and mod to extract all the digits of the number.
|
1 2 3 4 |
-- numTodigs :: Integral a => a -> [a] numTodigs 0 = [] numTodigs x = numTodigs (x `div` 10) ++ [x `mod` 10] |
the pad function is a simple recursive function that prepends ‘0’ until the length is as needed.
|
1 2 3 4 5 |
-- pad :: Integral a => [a] -> [a] pad xs | length xs >= 8 = xs | otherwise = pad (0 : xs) |
The mid4 function uses the substr function written in a previous post. i.e.
|
1 2 3 4 5 |
-- mid4 :: [a] -> [a] mid4 = substr 2 4 substr :: Int -> Int -> [a] -> [a] substr s e xs = take e (drop s xs) |
and finally digsToNum builds up the number, from a list, using powers of 10.
|
1 2 3 4 5 6 7 |
-- digsToNum :: [Int] -> Int digsToNum ds = eval (10 ^ l) ds where l = length ds - 1 eval _ [] = 0 eval n (x:xs) = (n * x) + eval (n `div` 10) xs |
So there it is, fairly trivial, but does highlight the notion of producing more complex functions by the composition of smaller, simpler functions.
To see it in action we need to use the iterate function to repeatedly apply vnRan like this:
|
1 2 3 4 5 |
-- λ-> take 100 $ iterate vnRan 1213 [1213,1369,4161,3921,4241,6081,8561,721,9841,5281,8961,9521,9441,2481,5361,321,3041,7681,7761,3121,641,881,6161,7921,2241,2081,561,4721,7841,1281,961,3521,7441,8481,7361,4321,1041,3681,9761,7121,8641,6881,8161,1921,241,8081,2561,8721,5841,7281,2961,7521,5441,4481,9361,8321,9041,9681,1761,1121,6641,2881,161,5921,8241,4081,4561,2721,3841,3281,4961,1521,3441,481,1361,2321,7041,5681,3761,5121,4641,8881,2161,9921,6241,81,6561,6721,1841,9281,6961,5521,1441,6481,3361,6321,5041,1681,5761,9121] *Main λ-> |
Here is the entire code in one block and I’ve added a main function that uses the system time to seed the vnRan with a ‘random’ 4 digit number.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
-- import Data.Time.Clock.POSIX numTodigs :: Integral a => a -> [a] numTodigs 0 = [] numTodigs x = numTodigs (x `div` 10) ++ [x `mod` 10] pad :: Integral a => [a] -> [a] pad xs | length xs >= 8 = xs | otherwise = pad (0 : xs) mid4 :: [a] -> [a] mid4 = substr 2 4 substr :: Int -> Int -> [a] -> [a] substr s e xs = take e (drop s xs) digsToNum :: [Int] -> Int digsToNum ds = eval (10 ^ l) ds where l = length ds - 1 eval _ [] = 0 eval n (x:xs) = (n * x) + eval (n `div` 10) xs sq :: Int -> Int sq = (^ 2) vnRan :: Int -> Int vnRan = digsToNum . mid4 . pad . numTodigs . sq sysInt :: IO Int sysInt = round . (* 1000) <$> getPOSIXTime fourDigits :: Int -> Int fourDigits n = n `mod` 10000 main = do n <- sysInt print $ take 150 $ iterate vnRan (fourDigits n) |
Anyone playing with this will soon see the sequence quite quickly either repeats or converges to zero and then repeats. This paper modifies the algorithm slightly and overcomes the quick convergence of the sequence and I’ll hopefully post about it soon!
Thanks for reading!