# Generating an Array of Random Floats

February 11 2007 by
false
Maple

2

While numerically testing a solution for the February 2007 IBM Ponder This challenge I had to generate a large n×2 Array of random floats uniformly distributed in [0,1).

My first attempt was to use structured flavors in the RandomTools package:

picks[4] := proc(n::posint)
uses RandomTools;
Array(Generate(list(list(float(range=0..1,'method=uniform'),2),n))
, datatype=float[8]);
end proc:

While this worked, it was rather slow. I then noticed that there is a listlist flavor in the RandomTools package. Using it gave a significant improvement (about 2×) over using nested list flavors.

picks[3] := proc(n::posint)
uses RandomTools;
Array(Generate(listlist(float(range=0..1,'method=uniform'),n,2))
, 'datatype=float[8]'
);
end proc:


By specifying a 2×n listlist structure and using the transpose option to the Array constructor this could be slightly improved:

picks[2] := proc(n::posint)
uses RandomTools;
Array(Generate(listlist(float(range=0..1,'method=uniform'),2,n))
, 'datatype=float[8]'
, 'transpose=true'
);
end proc:


While the structured flavors are convenient, they are not as fast as I'd like. My next attempt was to use the GenerateFloat64 procedure in the MersenneTwister subpackage of RandomTools. It calls (I believe) compiled code that is significantly faster than the random number generator used to generate structured flavors. Unlike the structured flavor approach, GenerateFloat64 returns only a singly random float that is (thankfully) uniformly distributed on [0,1) so we'll have to use an appropriate constructor to initialize the Array.

picks[1.3] := proc(n::posint)
uses RandomTools;
Array(1..n, 1..2, (i,j) -> MersenneTwister:-GenerateFloat64()
, 'datatype=float[8]'
);
end proc:
That is significantly faster (about 5× with n=50000) than the previous method. As you may have guessed by the numbering scheme, there are still a few improvements. Each is somewhat faster than the next.
picks[1.2] := proc(n::posint)
uses RandomTools;
local i;
Array([seq([MersenneTwister:-GenerateFloat64()
, MersenneTwister:-GenerateFloat64()]
,i=1..n)
]
, 'datatype=float[8]'
);
end proc:

picks[1.1] := proc(n::posint)
uses RandomTools;
local i;
Array([NULL
, [seq(MersenneTwister:-GenerateFloat64(), i=1..n)]
, [seq(MersenneTwister:-GenerateFloat64(), i=1..n)]
]
, 'transpose = true'
, 'datatype=float[8]'
);
end proc:


Jacques Carette suggests the following trick, which directly calls the compiled procedure used by GenerateFloat64.

picks[1.0] := proc(n::posint) local t;
uses RandomTools; kernelopts(opaquemodules=false):
t := (a,b) -> MersenneTwister:-MTKernelInterface(3);
Array(1..n, 1..2
, t
, 'datatype'=float[8]
);
end proc:

Here is a minor improvement. First, I restored the setting of opaquemodules. I also used the faster technique of passing a listlist structure rather than using functional initializer. Finally, to generate the lists quickly, two seq commands were used and the structure transposed, using the transpose option to Array.

picks[0.8] := proc(n::posint)
local i,rnd,opacity;
uses RandomTools;
opacity := kernelopts('opaquemodules'=false);
rnd := MersenneTwister:-MTKernelInterface;
kernelopts('opaquemodules'=opacity);
Array([NULL
, [seq(rnd(3), i=1..n)]
, [seq(rnd(3), i=1..n)]
]
, 'transpose=true'
, 'datatype=float[8]'
);
end proc:

This, too, can be further improved. For this case it is faster to create an empty Array and then use the map function to initialize random data.

picks[0.7] := proc(n::posint)
local opacity,rnd;
opacity := kernelopts('opaquemodules'=false);
rnd := RandomTools:-MersenneTwister:-MTKernelInterface;
kernelopts('opaquemodules'=opacity);
map(()->rnd(3), Array(1..n,1..2,'datatype=float[8]'));
end proc:

A bit more speed can be gained by using a hack. A number, including a hardware float, that is applied to any argument evaluates to itself. For example, 3.2(0.0) = 3.2. We can use this trick to avoid having to evaluate the operator ()->rnd(3), instead we use the unevaluated function 'rnd(3)' which, when evaluated, returns a hardware float that is then applied to 0.0 (the initial value of the Array) and so evaluates to itself:

picks[0.65] := proc(n::posint)
local opacity,rnd;
opacity := kernelopts('opaquemodules'=false);
rnd := RandomTools:-MersenneTwister:-MTKernelInterface;
kernelopts('opaquemodules'=opacity);
map('rnd(3)', Array(1..n,1..2,'datatype=float[8]'));
end proc:

Acer points out that rtable has a built-in frandom initializer. While the documentation for this feature is deficient (it doesn't mention whether the generated random numbers are uniformly distributed within the range, or other details) it is orders of magnitude faster than the previous approach.

picks[0.1] := proc(n::posint)
rtable(1..n, 1..2
, frandom(0..1, 1)
, 'subtype = Array'
, 'datatype=float[8]'
);
end proc:

Here is a performance comparison of all the techniques. and scaling the result. The fastest procedure is abpit 3 orders of magnitude faster than my first attempt. All tests were run with Maple 11 on a Linux box.

GGenerate a 50000 x 2 Array of random floats uniformly distributed in [0,1)

----- timing -------  ---------- memory ------------
total (s)  maple (s)  used (MB)  alloc (MB)  gctimes  proc
---------  ---------  ---------  ----------  -------  ----
0.02       0.02      0.76        0.81         0     picks[.1]
0.85       0.82     18.32       10.87         1     picks[.65]
1.03       0.92     18.96       10.69         1     picks[.7]
1.14       1.11     20.82       13.19         1     picks[.8]
1.25       1.20     21.97       14.62         1     picks[.85]
1.37       1.37     21.64       11.75         1     picks[.9]
1.50       1.46     21.74       11.75         1     picks[1.0]
1.41       1.37     23.43       14.31         1     picks[1.1]
1.53       1.50     24.12       16.06         1     picks[1.2]
1.68       1.61     23.41       13.25         1     picks[1.3]
7.79       7.58    125.29       39.74         4     picks[2]
8.51       8.36    129.21       41.49         4     picks[3]
16.63      16.01    199.61       43.93         6     picks[4]