To visually confirm these results, several images were created (see Figures 1-9). 256x256 pixel images were used to represent the first 64K values of the array. Each 32-bit value was broken into bytes, where byte 0 was used as the x-coordinate and byte 1 was used as the y-coordinate, with a white pixel drawn at this coordinate. The image was initialized to black.
The Windows rand() function shows a low density of values in Image 1 and patterns in Images 2 through 4, whereas Image 5 (15-bits per call) appears more random, visually confirming issues displayed in Table 1. Images 6 through 9 show ran2() which appears random for 1-bit, 2-bit, 3 and 4-bits extracted per function call.
The implementation in Listing 1 has a subtle flaw. Functions ran0(), ran1(), and ran2() generate a random number between 0.0 (inclusive) and 1.0 (exclusive) . This random number (0.0 to almost 1.0) is then multiplied by 65535, causing the value of 65535 (0xffff) to not be generated. Tables 1 through 4 show statistics for the following experiment. For 100 million random 16-bit values, a count of each possible value was gathered. Tables 1 through 4 show the minimum count value, the maximum count value, and the average count values. Not created column shows how many values had a count of zero; i.e., the random-number generator did not create that value. The left side of the tables shows the results for the flawed usage in Listing 1.
From the left side of Tables 1 through 3 several trends can be seen. When extracting 13-bits or more per call, ran0(), ran1(), and ran(2) do not generate evenly distributed all possible 16-bit values, with one value (0xffff) not generated at all by ran1() and ran2(). Table 4 shows an issue with the Windows C rand() function, where 6-bits or fewer do not create a significant number of possible values, thus skewing the statistics significantly.
To work around this problem, I used the following code:
float f1 = ranX( &randState ); // 0.0 .. 1.0 float f2 = ranX( &randState ); // 0.0 .. 1.0 unsigned long tmp1 = (unsigned long)( ( 1 << 16 ) * f1 ); unsigned long tmp2 = (unsigned long)( ( 1 << 16 ) * f2 ); unsigned long result = ( tmp1 << 16 ) | tmp2;
This fix is not necessary for the Windows C rand() function, since it generates 15-bits per call, and not floating-point numbers. Also, the ((1 << 16) - 1) multiplier is not used. This fix is also not necessary for Mersenne twister implementation, since it does not return a floating-point value.