Random Number Generation

1
Random Number Generation
2
Outline

• Problem statement
• Lehmers algorithm
• Period and full-period RNGs
• Modulus and multiplier selection (Lehmer)
• Implementation - overflow

3
Problem Statement

• We need to generate random numbers
• What does this mean?
• Successive draws are independent, unpredictable
• Each draw follows a certain probability
  distribution
• Here, interested in RNs uniformly distributed in
  interval (0, 1)
• Truly random numbers (whatever that means!) can
  be generated by hardware devices
• Measure time between successive decays of
  radioactive particles
• Measure white noise in electronic circuits
• Measure head movement on a rotating disk
• A problem with these (hardware) numbers is they
  truly are random, i.e., are not repeatable
• For many applications, pseudo-random number
  generators are used
• Nothing random about them only their output
  looks random

4
Desired Properties

• Randomness produce output satisfying statistical
  tests of randomness
• Repeatability ability to reproduce random number
  stream if necessary
• Portability easy to move random number generator
  (RNG) to a new machine
• Efficiency Efficient in time and memory
  requirements
• Documentation theoretical analysis and
  extensively tested

5
Random Number Generator

• Define a function
• u rand()
• that produces a floating point number u, where
  0ltult1 and any value greater than 0 and less than
  1 is equally likely to occur
• Exclude values 0.0 and 1.0 as possible outputs
  this will simplify things later when we compute
  functions based on u
• Restated
• For a large integer m, let ?m 1, 2, m-1
• Draw in integer x ? ?m at random
• Compute u x/m
• m should be very large
• Thus, our problem reduces to determining how to
  randomly select an integer in ?m

6
Lehmers Algorithm

• Generate a sequence of random integers in ?m
• x0, x1, x2, xi, xi1,
• Use the last drawn random integer to generate the
  next one xi1 g(xi) for some function g()
• g() is defined using two fixed parameters
• Modulus m, a large, fixed prime integer
• Multiplier a, a fixed integer in ?m
• and an initial seed x0 ? ?m
• Define g(x) a x mod m
• The mod function produces the remainder after
  dividing the first argument by the second
• More precisely u mod v u - v ?u/v? where ?x?
  (read floor) is the largest integer n such that
  n x

7
Observations

• xi1 a xi mod m
• The mod function ensures a value lt m is always
  produced
• Can the generator produce the value 0?
• If it does, all subsequent numbers in sequence
  will be zero
• It can be shown that if the m is prime, and the
  initial seed is non-zero, the generator will
  never produce the value 0
• Therefore this does produce values in ?m 1, 2,
  m-1
• The above simulates drawing balls from an urn
  without replacement, where each value in ?m
  represents a ball
• Actually, this violates our requirement for
  randomness that successive draws be independent
  (unpredictable)! We should be doing this with
  replacement
• Reasonable approximation if number of draws ltlt m
• Quality of the random number generator is
  dependent on good choices for a and m

8
Full Period Sequences

• Consider sequence produced by xi1 a xi mod m
• Once we repeat a value, the sequence repeats
  itself
• Sequence x0, x1, , xi, xip where xi xip
• p is called the period clearly p m-1
• In fact, it can be shown, that if we pick any
  initial seed x0, we are guaranteed this initial
  seed will reappear
• LP (Lemmis, Park, p. 42) theorem 2.1.2
• If x0 ? ?m and the sequence x0 x1 x2 is
  produced by the Lehmer generator xi1 a xi mod
  m where m is prime, then there is a positive
  integer p with p m-1 such that x0, x1, , xp-1
  are all different and
• xip xi i 0, 1, 2,
• In addition, (m-1) mod p 0
• Ideally, the generator cycles through all values
  in ?m to maximize the number of draws that are
  allowed, and guarantee any number can be produced
• Called a full-period sequence (p m-1)
• Non-full period sequences effectively partition
  ?m into disjoint sets

9
Modulus and Multiplier Selection

• Would like m to be as large as possible
• m 2i - 1 where i is the machine precision is
  the largest possible positive integer on a 2s
  complement machine
• Recall m must be prime
• It happens that 231-1 is prime
• Alas, 215-1 and 263-1 are not prime -(
• Would like full period sequencer (p m-1)
• For a given m, select multiplier a to achieve
  full period
• Algorithm to test if a is a full period
  multiplier (m must be prime)

p 1 x a // assume, initial seed is
1 while (x ! 1) // cycle through numbers until
repeat p x (a x) m // careful
overflow possible if (p m-1) // a is a full
period multiplier else // a is not a full
period multiplier
10
Other Useful Properties

• Theorem 2.1.1LP If the sequence x0, x1, x2,
  is produced by a Lehmer generator with multiplier
  a and modulus m, then
• xi ai x0 mod m, i0, 1, 2,
• Note this is not a good way to compute xi!
• Theorem 2.1.4LP, p. 45 If a is any full-period
  multiplier relative to the prime modulus m, then
  each of the integers
• ai mod m ? ?m i1, 2, 3, ... m-1
• is also a full period multiplier relative to m if
  and only if the integer i has no prime factors in
  common with the prime factors of m-1 (i.e., i and
  m-1 are relatively prime)

// Given prime modulus m and any full period
multiplier a, // generate all full period
multipliers relative to m i 1 x a //
assume, initial seed is 1 while (x ! 1) //
cycle through numbers until repeat if
(gcd(i,m-1)1) // xaimod m is full period
multiplier i x (a x) m // careful
overflow possible
11
Implementation

• Assume m231-1
• Problem
• Must compute a x mod m
• Obvious computation is compute (a x) first, then
  do mod operation
• The multiplication might overflow, especially if
  m-1 is large!
• Floating point solution
• Could do arithmetic in double precision floating
  point if multiplier is small enough
• Double has 53-bits precision in IEEE floating
  point
• May have trouble porting to other machines
• Integer arithmetic faster than floating point

12
Another Solution

• Compute a x mod m
• Can we do mod operation first, before multiply?
• Basic idea suppose it were the case that m a q
• a x mod m a x mod a q a (x mod q)
• x mod q is at most q-1
• Thus a (x mod q) at most a(q-1) lt aq m no
  overflow!
• Of course, m is prime, so let m a q r
• q quotient r remainder
• Let
• ?(x) a (x mod q) - r ?x/q?
• ?(x) ?x/q? - ?ax/m?
• It can be shown
• ax mod m ?(x) m ?(x)
• and ?(x) 0 if ?(x) ? ?m , and ?(x) 1 if -?(x) ?
  ?m
• Algorithm maqr is prime, rltq, and x? ?m
• t a (x q) - r (x / q) // ?(x)
• if t gt 0 return t
• else return (t m)

13
Concluding Remarks

• Random number generation a rich field still an
  active research area
• No consensus on the best random number
  generator (but Lehmers algorithm with
  appropriate choice of a and m is considered a
  good one)
• Many RNG have been proposed in the literature
• Lehmers algorithm an example of linear
  congruential type of generator
• Pick multiplier and modulus based on principles
  (like full-period) and do extensive statistical
  tests to validate randomness
• Other issues arise, e.g., if a very long sequence
  is needed (e.g., long running simulations)