Streams of Pseudo-Random Numbers in the Functional Programming Language Haskell 1.1 2006/04/02 15:12:55.370 GMT-5 2006/04/09 22:23:21.374 GMT-5 Jacob Schrum he_deceives@hotmail.com Jacob Schrum he_deceives@hotmail.com functional haskall lists numbers programming pseudo random streams 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) whereThe 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. Inputs 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) Inputs 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. Inputs 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` rangeThe 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. Inputs 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. Inputs 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. Inputs minimum value maximum value seed >boundRandomFloatings :: Floating a => a -> a -> Integer -> [a] >boundRandomFloatings lo hi n = map ((+ lo).(* (hi - lo))) (smallRandomFloatings n)