In the past, I’ve tended to use the Mersenne Twister for generating pseudo-random numbers. This is a widely used generator because it has a very long period of 2^{19937} – 1. The one issue is that it’s not the fastest generator and this is particularly an issue when drawing random numbers for implementing the Gillespie direct method. A much faster variant is the Multiply-with-carry (MWC) random number generator. This uses much simpler arithmetic which means it’s faster. The algorithm can be described in three lines:

The symbol >> means shift the bits to the right by the given number and << is similar but moves the bits to the left. In object pascal, these operations are represented by shr and shl respectively. The generator needs two values to be initialized, m_z and m_w I've recently implemented the Gillespie method in three languages, C, Go and Object Pascal. Object Pascal has good support for the Mersenne Twister but I had to create my own MWC generator. The code below shows the object pascal implementation:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 |
// ------------------------------------------------------------ // George Marsaglia's MWC (multiply with carry) generator. // This code translated from C# to Object Pascal and obtained from // https://www.codeproject.com/Articles/25172/Simple-Random-Number-Generation unit uRandMWC; interface Uses System.Types; // Call this before doing anything procedure initMWC; // If you want to change the seed call this procedure setSeed (seed : cardinal); // Call this to get a random number between 0 and 2^32 function getUint() : cardinal; // Call this to get a uniformly distributed random number between 0 and 1.0 function getUdouble : double; implementation // In Delphi 32 and 64 bit, cardinal is an unsigned 32-bit integer var m_z, m_w : cardinal; procedure setSeed (seed : cardinal); begin m_z := seed; m_w := 678934; end; procedure initMWC; begin m_z := 467567; m_w := 125681; end; function getUint() : cardinal; begin m_z := 36969 * (m_z and 65535) + (m_z shr 16); m_w := 18000 * (m_w and 65535) + (m_w shr 16); result := (m_z shl 16) + m_w; end; function getUdouble : double; var u : dword; begin u := getUint(); // The magic number below is 1/(2^32 + 2). // The result is strictly between 0 and 1. result := double (u + 1.0) * 2.328306435454494e-10; end; end. |

What about some tests? There are random generator test systems such as TestU01 and dieharder but they are written to run on Linux and I didn't have time to set something up on windows plus the usage docs are not that great. Instead I did three simple sanity tests. The simplest is to check that the mean value for a sample of random numbers is 0.5. Another test is to make sure that, say out of 10,000 random numbers, 10% are between 0 and 0.1, etc. Ideally there should be statistical tests done to check that these values are within certain tolerances. The third test is a visual test where I generate 32768 random numbers and map them on to a bitmap. Each random number is 32 bits long and each bit represents a pixel in the bitmap. 32768 numbers means I generate a bitmap 1024 by 1024 pixels. Thus 32 numbers occupy each row of the bitmap. I wrote a simple Delphi console application that calls the random number generator and creates the bitmap. The main part of this code is shown below. For a comparison I also implemented the NOT to be used randu generator and also the Delphi's built in generator. The approach I take isn't very sophisticated and I am sure it could be made shorter, but for each random number I generate, I convert it into a string of 1's and 0's. I then work my way through the string setting pixels in a bitmap. I do this 32 times for each row in the bitmap.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 |
program RandomGen; {$APPTYPE CONSOLE} {$R *.res} uses System.SysUtils, System.Types, Vcl.Graphics, uRandMWC in 'uRandMWC.pas'; var r1 : double; ri : cardinal; i, j, k : integer; bmp : TBitmap; sb : string; iCol : integer; state : cardinal; count : integer; // Get a single bit as position Bit from cardinal aValue function getBit(const aValue: cardinal; const Bit: Byte): Boolean; begin Result := (aValue and (1 shl Bit)) <> 0; end; // Converts a cardinal (32-bits) into a string of ones and zeros function getBinaryString (aValue : cardinal) : string; var i : integer; begin result := ''; for i := 0 to 32 - 1 do if getBit (aValue, i) then result := result + '1' else result := result + '0' end; function randu : cardinal; begin state := (65539 * state) mod 2147483648; result := state; end; begin state := 25; // starting value for radu writeln ('wait.....'); initMWC(); bmp := TBitmap.Create; bmp.PixelFormat := pf1bit; bmp.Width := 1024; bmp.Height := 1024; count := 0; for i := 0 to bmp.Height - 1 do begin iCol := 0; for j := 0 to bmp.Width div 32 - 1 do begin // Choose a random number generator //ri := randu; ri := getUint(); //ri := random (4294967296-1); sb := getBinaryString (ri); for k := 0 to 32 - 1 do begin if sb[k+1] = '0' then begin bmp.Canvas.Pixels[iCol,i] := clWhite; end else begin bmp.Canvas.Pixels[iCol,i] := clBlack; end; inc (iCol); end; end; inc (count, bmp.Width div 32); end; bmp.SaveToFile('image.bmp'); writeln ('count = ', count); writeln ('Press any key to coninue'); readln; end. |

The results of running this code is below. First randu to show you how bad it is, note the patterns and lines in the bitmap:

The next figure shows the results from Delphi’s builtin random number generator, which isn’t so great either:

Finally, the MWC algorithm, which you can see from the image the distribution of pixels appears to show no pattern: