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 ‘

**‘ then ‘**

*convert number to digits***‘ then ‘**

*pad to 8*

*get***‘ then ‘**

*middle 4***‘**

*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!