Þ   briarpig  » thorn  » demos  » rand

     Pseudo random number generators are useful for several things, and Wil often uses one standard implementation based on Park and Miller's Oct1988 Random number generators: good ones are hard to find, in which they identify three different algorithms with good statistical behavior when implemented in a way avoiding 32-bit overflow despite 32-bit arithmetic.

     One of the nice properties of this generator is the period: it visits every integer value from one to 2*31-2 exactly one time before repeating. Among other things, this can be used to generate floating point random numbers with a uniform distribution between zero and one without ever being equal to zero or one.

     Variant r64rand() below is nearly identical, differing only in converion to double and division by 0x7FFFFFFF which is one more than the largest value returned by h32rand().

     As a result, r64rand() returns values strictly more than zero and strictly less than one, distributed uniformly between zero and one (because the numerator uses all integers between zero and 0x7FFFFFFF: each equally likely).

     Note neither r64rand() nor h32rand() is the fastest way to calculate these answers. Wil just decided to stick with the variant shown in Park and Miller's paper as long as performance was not the greatest concern. Sometimes Wil uses h32rand() for hashing, but most of the time Wil only uses these methods in monte carlo testing with pseudo random data sets.

     Below is the original version of code for Everett F. Carter's Box-Muller transformation. The yr64g class (top column right) is based directly on this code. This was originally posted on treedragon in March 2001, but you can find a copy on this site in Aug2007. See yr64g for an explanation.

     See http://www.taygeta.com/random/gaussian.html for Dr. Carter's explanation. But the main idea you should get is this: box_muller() returns values distributed in a normal distribution in the shape of a bell curve. The mean and standard deviation are based on args you pass to the function.

     Dr. Carter's conditions for use are specified above.

     Most code on this site is available only under the BriarPig mu-babel license described fully on the rights page. But the random number generators on this page are an exception: if you remove the y prefix and attribute "BriarPig was here", you can use the BSD license as a dual license, for code on this page only.

     You do not have permission to reprint this page in any way. No feeds or repackaging allowed. You can link this page if you want folks to read it. (RSS feeds are forbidden, in case you have trouble working out that special case.)

     Sample code below shows use of method ydvp::pjigsaw() listed bottom right, one column over, which constructs pseudo random partitions of byte sequences in yv format (cf «).

     Class ydvp appears in the iovec demo (cf «) representing an iovec array in iterator form; yd256vp is a variant with 256 iovecs built-in (cf «). The following code uses a yd256vp instance to hold a randomly constructed iovec array partitioning a string into contiguous, adjacent, disjoint pieces.

     Each time ydvp::pjigsaw() is called, it creates a one-dimensional jigsaw puzzle of vector lohi which still contains all the same bytes of lohi, and in the same order, but partitioned into a different set of iovecs whose sizes are pseudo randomly assigned with a target average size and standard deviation.

     (Why would anyone want to do this? An example appears in a future compare demo, where Wil tests code comparing content in one iovec array to another in a generalized scatter gather pattern finding method seeking long matches. Known but scrambled matches can be generated with ydvp::pjigsaw() then searched to see if expected matches are found.)

     Here's the output of code above:

     Once Wil gets a unit test working well with small sample sets, he usually cranks up the number of iterations to more than a hundred thousand test cases, to ensure many variations of odd size combinations are exercised. Wil doesn't want to find out later that some edge case had a bug in single byte iovecs (for example) in certain positions when multiple coincidences occur.

     With pseudo random monte carlo testing, you might as well run an obscenely large number of sample data set cases if the elapsed time isn't enough to keep a developer waiting. (And even then, the developer can wait anyway for i/o tests that can't be rigorously covered in less than a few minutes of wall clock time.) In-memory unit tests are over in a flash.

     thorn: todo, names, fd, iovec, assert, log, run, hex, crc, buf, in, out, quote, escape, compare, file, deck, cow, arc, blob, tree, slice, rand « Þ, time, stat, hash, heap, node, primes, page, book, pile, stack, atomic, lock, mutex, thread, map, meter, list, iter, ctype

     (mu: toy, peg, imm, tag, box, symbol, token, number, bigint, class, method, reader, writer, eval, env, vm, gc, world, pcode, compiler, asm, lathe, lisp, smalltalk, design, weight, jar, card, harp, debug, profile)

     Some demos are stubs: todo is a demo guide. See toy for mu updates on language pages; names introduces naming schemes.

     Class yr64g below is Wil's fourth or fifth rewrite of Carter's box_muller() (bottom column left). Name yr64g means thorn real 64-bit generator. The main function grandom() returns the same value as Carter's box_muller(). Method gaussian() is an inline synonym since normal distribution is also called gaussian distribution; operator++() is yet another synonym.

     An interesting demonstration of gaussian generator yr64g appears in the future stat demo, which empirically verifies the box-muller tranformation really does emit a distribution with the requested mean and standard deviation. (The yrsum class in the stat demo calculates the mean and standard deviation in a sequence of values.)

     Before diving into the code, let's consider a couple of your potential questions. How does the box-muller tranformation work? And why would you want gaussian bell curve distributions in random number generators? What's it good for?

     First, Wil doesn't know how it works; you can read the literature if you really want to know. Wil just treats the algorithm as a little bit of convenient magic simplifying test data set generation. After empirically seeing it do the right thing in the stat demo, Wil no longer has concerns about how it works. It just does.

     Second, bell curve distributions are great to model natural data sets of variable sizes that still cluster around expected values with some typical variance from the average value. Wil first used a random gaussion distribution to model what happens when caching web pages in a reverse proxy server, given a population of web pages having an expected average size and consistent standard deviation.

     Among other things, this is a good way to exercise algorithms that optimize data sets with locality better than totally random data sets. The pattern of clustering around the average at the top of a bell curve is a good appoximation of data access patterns with some locality.

     A seed to constructors passes by pointer so yr64g can share an external random number stream with other users.

     This version of grandom() closely follows Carter's original.

     The remaining two methods are applications of the Park and Miller random number generators; guniform() is central to grandom() and ought to be inlined (if speed matters).

     Code for ydvp::pjigsaw() almost appeared earlier in the iovec demo where ydvp was defined (cf «). But since that demo was already very long, this method using yr64g was saved for presentation here.

     Sample code showing typical effect for an input yv sequence is shown bottom column left (cf «). So let's just start with a documentation comment in the code explaining what it does. (Wil often explains complex methods with comments when intended effect is otherwise hard to grasp, as in this case.)

     "Dang it, Wil," complained Stu. "Can't you write a slightly more crisp summary of that here?"

     "Yes, I could," replied Wil. "But you're getting spoiled. That description's not half bad. Did you try to read it?"

     "I think I have an acquired aversion to reading code comments," Stu explained. "I'm gunshy. Have pity on me."

     "Oh you spazz," Wil sighed. "Start at Box-Muller in the comment and read each sentence slowly. You'll get it."

     This opening part of the method first checks whether the iovec array length p_x is zero, and does nothing given no iovecs. Then it creates a copy of the random seed in local seed in case some other thread is also changing it. Then it checks whether the source yv vector is too small to bother partitioning, or whether the iovec array is too small to aim to use only half; either way, all of src becomes one iovec as output.

     Here pjigsaw() aims for an average fragment size that would use half the iovecs in the ydvp iovec array, with a standard deviation of two thirds the average. But eight is the min for average.

     These pointers let pjigsaw() walk both source yv and destination iovec array while reserving the last iovec as a safety catch-all iovec, in case random chance uses all the iovecs with tiny fragments.

     The rest should be clear enough from earlier comments. Final iovec array length is set on the line assigning p_n toward the end.