Reversible Random Number Generator

November 24, 2015

Since we have a modular inverse function from a previous exercise, this is easy:

(define next #f)
(define prev #f)
(let ((a 69069) (c 1234567) (m (expt 2 32)))
  (set! next (lambda (x) (modulo (+ (* a x) c) m)))
  (set! prev (lambda (x) (modulo (* (inverse a m) (- x c)) m))))

The two defines that are re-bound inside a let are the standard Scheme idiom for keeping two variables private to a cooperating set of function. Here are some examples:

> (next 2718281828)
3103402651
> (next 3103402651)
4281062310
> (next 4281062310)
1670430837
> (prev 1670430837)
4281062310
> (prev 4281062310)
3103402651
> (prev 3103402651)
2718281828

You can run the program at http://ideone.com/mRuFZS, where you will also see the function that performs the modular inverse. If you don’t want to use a linear congruential generator, an alternate method of solving the problem is to encrypt the numbers 1, 2, 3, … using your favorite encryption method to convert the sequential index to a random number.

Advertisements

Pages: 1 2

2 Responses to “Reversible Random Number Generator”

  1. Globules said

    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 123456
    
    $ ./revrand 
    mmix: success!
    vms: success!
    vb: success!
    minstd_rand: success!
    bogus: 1234 and 123456 are not coprime.
    
  2. Jussi Piitulainen said

    This 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
    
    end
    

    Julia 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)
    

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 )

Google photo

You are commenting using your Google 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 )

Connecting to %s

%d bloggers like this: