Reversible Random Number Generator
November 24, 2015
Some applications of random number generators. games, for instance, require that the random sequence be available to run in reverse.
This is easy to do with a simple linear congruential random number generator, which is characterized by the formula next = a * prev + c (mod m). With a little bit of algebra, that formula can be written prev = a−1 * (next - c) (mod m), where a−1 is the inverse of a (mod m); that modular inverse always exists because the linear congruential generator requires that a and m are coprime.
Your task is to write a reversible random number generator. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
A Haskell solution.
import System.Random import Text.Printf -- Extended Euclidean algorithm. Given non-negative a and b, return x, y and g -- such that ax + by = g, where g = gcd(a,b). Note that x or y may be negative. gcdExt :: Integral a => a -> a -> (a, a, a) gcdExt a 0 = (1, 0, a) gcdExt a b = let (q, r) = a `quotRem` b (s, t, g) = gcdExt b r in (t, s - q*t, g) -- Given a and m, return Just x such that ax = 1 mod m. If there is no such x -- return Nothing. modInv :: Integral a => a -> a -> Maybe a modInv a m = let (i, _, g) = gcdExt a m in if g == 1 then Just (pos i) else Nothing where pos x = if x < 0 then x + m else x type LCG a = a -> a -- Given a, c and m return Just (fwd, rev), where fwd is a forward linear -- congruential generator and rev is its corresponding reverse generator. -- Return Nothing if a and m are not coprime. lcgPair :: Integral a => a -> a -> a -> Maybe (LCG a, LCG a) lcgPair a c m = case a `modInv` m of Nothing -> Nothing Just a' -> let fwd prev = (a*prev + c) `mod` m rev next = a'*(next - c) `mod` m in Just (fwd, rev) -- Test the LCG pairs by generating one sequence of numbers using the forward -- generator, then a second sequence using the reverse generator, starting -- with the last number produced by the forward generator. Compare the first -- sequence with the reverse of the second one, indicating success if they are -- the same, or failure if they are not. testLcg :: String -> Integer -> Integer -> Integer -> Integer -> IO () testLcg name seed a c m = case lcgPair a c m of Nothing -> printf "%s: %d and %d are not coprime.\n" name a m Just (fwd, rev) -> let fs = take 10000 $ iterate fwd seed rs = take 10000 $ iterate rev (last fs) in if fs == reverse rs then printf "%s: success!\n" name else printf "%s: failure!\n" name main :: IO () main = do -- Ensure the seed is less than the smallest (valid) modulus below. seed <- randomRIO (0, 2^24-1) -- Test a few LCGs listed on the Linear congruential generator page -- (https://en.wikipedia.org/wiki/Linear_congruential_generator). testLcg "mmix" seed 6364136223846793005 1442695040888963407 (2^64) testLcg "vms" seed 69069 1 (2^32) testLcg "vb" seed 1140671485 12820163 (2^24) testLcg "minstd_rand" seed 48271 0 (2^31-1) testLcg "bogus" seed 1234 0 123456This implements the Julia (v0.4.1) iterator protocol (methods for start, done, next). I believe each version of the generator (Gen, ForwardIter, ReverseIter) is actually stored as a single 64-bit value as long as the compiler can infer its type. So the whole thing could even be held in a register. A call to next returns both the generated value and the next state. There are no stateful objects in this implementation.
module ReversibleGenerator export Gen, ForwardIter, ReverseIter import Base: start, done, next immutable Gen # 64-bit state n::UInt64 end immutable ForwardIter gen::Gen end immutable ReverseIter gen::Gen end # Linear congruential generator parameters # of "POSIX [jm]rand48" from Wikipedia, # reverse formula from Praxis. # m = 2^48, c = 11, return bits 47..15 (a 32-bit integer) a = 0x5deece66d b = 0xffffdfe05bcb1365 # = invmod(a, m) start(i::ForwardIter) = i.gen done(::ForwardIter, ::Gen) = false next(::ForwardIter, x::Gen) = let (UInt((x.n >>> 15) & 0xFFFFFFFF), # 32-bit value as system UInt Gen(mod(a * x.n + 11, 2^48))) end start(i::ReverseIter) = i.gen done(::ReverseIter, ::Gen) = false next(::ReverseIter, x::Gen) = let (UInt((x.n >>> 15) & 0xFFFFFFFF), # 32-bit value as system UInt Gen(mod(b * (x.n - 11), 2^48))) end endJulia for-loops use the iterator protocol implicitly, but then one gets access only to the values. To reverse the process in the middle requires access to the state, hence the explicit calls to start and next (but never done).
using ReversibleGenerator xs = Array(UInt, 0) # empty vector of UInt (system size) # i serves as a type for start/done/next # s is an actual generator/state i = ForwardIter(Gen(31415926)) s = start(i) for k in 1:3 x, s = next(i, s) push!(xs, x) end i = ReverseIter(s) for k in 1:4 x, s = next(i, s) push!(xs, x) end # prints a palindromic random array (in 0-padded hex :) println(xs)