**Algorithm 3.6:** Use the linear congruential method to generate a uniform set of pseudo-random numbers.

This is great and something I never did before – like most of the stuff here. The method were using is the *linear congruential method* which should create a uniform distribution of numbers.

The basic idea is to use the following formula to create the next pseudo-random number:

for n >= 0

However, there are specific criteria for each parameter.

*Parameter *:

*Parameter m*: m >= length sequence required, (ax+b) mod m is an integer

*Parameter a*: if m is a power of 2 then a mod 8 = 5; if m is a power of 10 then a mod 200 = 21

sqrt(m) < a < m – sqrt(m)

(a-1) should be a multiple of every prime dividing into m

if m is a multiple of 4 then (a-1) should be a multiple of 4

*Parameter b*: b should be odd, not a multiple of 5

These are a lot of requirements and the easiest thing to do is use existing values. Wikipedia provides common parameters.

The ones provided by the book are:

m = 4096

b = 853

a = 109

Here’s the code:

1 2 3 4 5 6 7 8 9 10 11 |
(defn gen-pseudo-random " Generates a sequence of random integers less than 4096." [seed] {:pre [((complement neg?) seed), (< seed 4096)]} (let [a 109 b 853 m 4096 random-number (mod (+ (* a seed) b) m)] (lazy-seq (cons random-number (gen-pseudo-random random-number))))) |

Super straight forward. Let’s see how good it works. We want to compare the average.

**Problem 3.6.1:** Confirm that the algorithm repeats after generating m random numbers. Compute the mean value and variance for the set of m pseudo-numbers.

1 2 3 4 |
user=> (last (take 4096 (gen-pseudo-random 40))) 40 user=> (last (take 4096 (gen-pseudo-random 2345))) 2345 |

Seems to work.

1 2 3 4 5 6 |
user=> (def pseudo-set (take 4096 (gen-pseudo-random 222))) #'user/pseudo-set user=> (average pseudo-set) 4095/2 user=> (variance pseudo-set) 1398101.25 |

Looks fine. Let’s see how big the variance over the average over the set is.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
user=> (dotimes [n 20] (println "Seed:" n "Average:" (average (take 4096 (gen-pseudo-random n))))) Seed: 0 Average: 4095/2 Seed: 1 Average: 4095/2 Seed: 2 Average: 4095/2 Seed: 3 Average: 4095/2 Seed: 4 Average: 4095/2 Seed: 5 Average: 4095/2 Seed: 6 Average: 4095/2 Seed: 7 Average: 4095/2 Seed: 8 Average: 4095/2 Seed: 9 Average: 4095/2 Seed: 10 Average: 4095/2 Seed: 11 Average: 4095/2 Seed: 12 Average: 4095/2 Seed: 13 Average: 4095/2 Seed: 14 Average: 4095/2 Seed: 15 Average: 4095/2 Seed: 16 Average: 4095/2 Seed: 17 Average: 4095/2 Seed: 18 Average: 4095/2 Seed: 19 Average: 4095/2 |

Yup. By the way, here are the functions for calculating the average and variance:

1 2 3 4 5 6 7 8 9 10 11 12 13 |
(defn average "Calculates the average of a collection" [coll] (/ (reduce + coll) (count coll))) (defn variance "Calculates the variance of a collection" [coll] (* (/ 1 (count coll)) (let [my-average (average coll)] (reduce + (map #(Math/pow (- % my-average) 2) coll))))) |

Pretty standard.