Analogmachine Blog: Marsaglia’s MWC (multiply with carry) Random Number Generator for Object Pascal

Originally Posted on December 14, 2019 by hsauro

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:

$m_z = 36969 * (m_z \text{ and } 65535) + (m_z >> 16)$

$m_w = 18000 * (m_w \text{ and } 65535) + (m_w >> 16)$

$\text{random number} = (m_z << 16) + m_w$

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.