Summary: In developing an Artificial Life simulation, a pseudo-random number generator for the functional programming language Haskell is designed with complete code provided.

Random numbers can be quite a controversial subject. One can argue that it is possible to predict everything given enough knowledge. That would mean that randomness is simply an abstraction for things we cannot understand, or for which there is not enough time or processing power to make carrying out the necessary calculations worthwhile, even if possible.

In computers the issue of randomness is somewhat less philosophical but still important. Strictly speaking, computers are completely deterministic. This means that generating random numbers is not possible. However, it is possible to generate pseudo-random numbers, which share many of the properties of random numbers and for nearly all applications are equally useful.

In most imperative languages programmers have the option of creating their own random number generators with seeds, which are initial numbers that influence the series of random numbers created by the generator. Some languages simplify this further and use the value from the computer's internal clock as a seed. The problem with this practice is that one cannot replicate an experiment if the clock is used for seeding.

In Haskell the deterministic nature of computers becomes even clearer. Haskell also has a built in random number generator, but because it is dependent upon the state of the computer it returns an IO type. This makes computations requiring a large number of random numbers difficult without the use of `unsafePerformIO`

, which converts an IO type to a pure type. In the interest of maintaining the functional integrity of code, this function is not used. Instead, a custom random number generator taken from the book "The Craft of Functional Programming" by Simon Thompson (1999) is used.

The method presented in Thompson's book is called the "linear congruential method," and all of the random number streams presented in this module stem from this method.

```
>module RandomStream
> (nextRandomIntegral,
> randomIntegrals, boundRandomIntegrals,
> smallRandomFloatings, boundRandomFloatings) where
```

The three constants below come directly from Thompson's book, and each plays a role in the generation of pseudo-random numbers. The `nextRandomIntegral`

function below takes an `Integral`

value as input and multiplies it by the constant named `multiplier`

. This value is then added to the `incrementer`

, and modulus division by the constant `modulus`

is performed on the result. The value `modulus`

is important because it sets the range for the function. All numbers generated by this function have a range [0, `modulus`

). The `modulus`

is set to 65536, which is more than large enough for all applications to which this random number generator is applied.

```
>multiplier :: Num a => a
>multiplier = 25173
>incrementer :: Num a => a
>incrementer = 13849
>modulus :: Num a => a
>modulus = 65536
```

According to Donald E. Knuth's book "The Art of Computer Programming Vol. 2: Seminumerical Algorithms," (1969) Thompson's choice of "magic numbers" is valid. Knuth's book deals with, among other things, the issue of pseudo-random number generation by computers. Knuth explains that the linear congruential method always creates a looping sequence of numbers, and while it would be desirable to have an infinite list of random numbers, one must instead settle for creating a sequence with as large a period as possible. The period is the number of values generated before looping occurs.

The maximum possible value of the period is `modulus`

, in which case every value in the sequence occurs exactly once before repeating. In order for the period to equal `modulus`

, several conditions must be met, as according to the Theorem A in Knuth's book:

`incrementer`

must be relatively prime to`modulus`

.- (
`multiplier`

- 1) is a multiple of p, for every prime p dividing`modulus`

- (
`multiplier`

- 1) is a multiple of 4, if`modulus`

is a multiple of 4.

Thompson's numbers satisfy condition 1 because gcd(13849,65536) = 1. Condition 3 is satisfied because 4|(25173 - 1) and 4|65536. Proving condition 2 is more difficult, but thankfully Knuth provides a simpler test that can be used to assure that `multiplier`

and `modulus`

are chosen so that every value in the range [0,`modulus`

) occurs before repeating. Since 65536 is a power of 2 (65536 = 2^16), all that must be done to assure that every number in the sequence occurs before repeating is to choose `multiplier`

such that `multiplier`

mod 8 = 5. Thompson's `incrementer`

satisfies this property (25173 mod 8 = 5).

One advantage of having the period equal `modulus`

is that it does not matter what value is used to seed the generator. Since all values from 0 to (`modulus`

- 1) appear in the sequence exactly once, any of these values can be used as a starting point. This is convenient because multiple copies of this generator are meant to be used simultaneously, so each needs to have a different starting seed in order to be independent from the others.

Now on to the actual random number generator. Notice that the type signature of `nextRandomIntegral`

allows for the generation of any general `Integral`

type. This means the function works for both `Int`

and `Integer`

types, which is convenient since random numbers of both types are needed by the simulation.

- seed or previous random number

```
>nextRandomIntegral :: Integral a => a -> a
>nextRandomIntegral n = ((n * multiplier) + incrementer) `mod` modulus
```

The `nextRandomIntegral`

function takes a single value as input and generates another. In order to obtain many random numbers this function needs to be applied repeatedly, and Haskell happens to have a useful function named `iterate`

which accomplishes this. Iterating with the `nextRandomIntegral`

function means that the result of each application of `nextRandomIntegral`

is itself used as the input to the next application of the function. The outputs are returned in a stream. This is possible because of Haskell's use of lazy evaluation. All that is needed is an initial value to start with, which is the seed of the random number generator. The seed is the input to the function `randomIntegrals`

and different seeds create different random streams. According to Thompson the numbers in this sequence "all occur with the same frequency." (Thompson 1999, pg. 367)

- seed

```
>randomIntegrals :: Integral a => a -> [a]
>randomIntegrals = iterate nextRandomIntegral
```

Thompson also provides a method of restricting the range of this function to generate numbers in the range [s,t]. This is done with the `scaleSequence`

function. It works by splitting the original range of numbers into "range blocks" and then mapping values in given ranges to new values. Unfortunately, it only works properly if the size of the new range is divisible by `modulus`

. Otherwise it maps the last few values in the initial range [0, `modulus`

- 1] to values outside of the desired range [s,t]. Therefore the function `boundRandomIntegrals`

removes values out of range after applying `scaleSequence`

to the initial stream of random numbers produced by `randomIntegrals`

.

- minimum value
- maximum value
- list of random integrals in range [0,
`modulus`

- 1]

```
>scaleSequence :: Integral a => a -> a -> [a] -> [a]
>scaleSequence s t = map scale
> where
> scale n = n `div` denom + s
> range = t - s + 1
> denom = modulus `div` range
```

The `boundRandomIntegrals`

function combines the `randomIntegrals`

and `scaleSequence`

functions into one simple function, and then removes the occasional value that is outside of the scaled range. All values within the desired range are equally likely to occur even though use of the `scaleSequence`

function may result in some values being discarded. The probability of `scaleSequence`

returning a value that needs to be discarded is slightly less than the probability of any particular valid value being returned. After invalid values are discarded the resulting distribution of possible random numbers is uniform. Calling this function with a desired range [a,b] generates a stream of numbers in this range. The seed for the random number generator is also needed.

- minimum value
- maximum value
- seed

```
>boundRandomIntegrals :: Integral a => a -> a -> a -> [a]
>boundRandomIntegrals a b n = filter (<= b) (scaleSequence a b (randomIntegrals n))
```

Random floating point numbers are also needed in the simulation. Thompson does not provide a method for generating random floating point numbers, but it is quite easy to derive a floating point number generator from the integral number generator he provides. First a function named `smallRandomFloatings`

that generates floating point numbers in the range [0,1] is defined. This generator is particularly useful on its own, and also serves as the basis for the function `boundRandomFloatings`

further below. The `smallRandomFloatings`

function works by dividing every number in the `randomIntegrals`

sequence by the maximum possible value in that sequence, namely `modulus`

- 1. This allows for `modulus`

number of possible random numbers in the range [0,1], which is more than enough for the purposes of the simulation.

- seed

```
>smallRandomFloatings :: Floating a => Integer -> [a]
>smallRandomFloatings n = map ((/ (modulus - 1)).fromInteger) (randomIntegrals n)
```

Since random floating point numbers larger than 1 and less than 0 are also needed in the simulation we need a function that scales floating point numbers to this range. Unlike the `boundRandomIntegrals`

function which scales the range of `randomIntegrals`

down, the `boundRandomFloatings`

function below scales the range of `smallRandomFloatings`

up. It can map values from this function to any range, but it is important to remember that no matter what the range, there are still only `modulus`

number of possible values. This means that the number generator below works best for small ranges in which high precision is not important. As such it is ideal for this simulation. However, this generator would not work well for generating floating point numbers in particularly large ranges, as increasing the range increases the space between consecutive values that this function is capable of generating. In fact, if the range is set to [0, `modulus`

- 1] then the space between each possible value is exactly one. In other words, the usefulness of the function `boundRandomFloatings`

decreases as the range increases. The function works by multiplying each value of the `smallRandomFloatings`

stream by the range, thus creating a stream with values from 0 to a - b, where a is the minimum and b is the maximum. The values of this stream are then shifted by being added to the minimum value. This makes the smallest possible value equal to the minimum of the range, even if it was negative, and the greatest possible value equal to the maximum of the range.

- minimum value
- maximum value
- seed

```
>boundRandomFloatings :: Floating a => a -> a -> Integer -> [a]
>boundRandomFloatings lo hi n = map ((+ lo).(* (hi - lo))) (smallRandomFloatings n)
```