Random numbers provided by "function Random(l:int64):int64;" are not equidistributed.
Original Reporter info from Mantis: Pangea
-
Reporter name: Pangea
Original Reporter info from Mantis: Pangea
- Reporter name: Pangea
Description:
The output of "function Random(l:int64):int64;" in system.inc is not equidistributed. The first 2^64 mod l numbers have a slightly larger probability to be chosen than the rest of the numbers. This is most noticeable when looking at probabilities like p(r < l/2) where r = Random(l) and l is large.
Steps to reproduce:
Consider the following program (expected output 1/2; actual output 2/3):
program random_test;
{$mode objfpc}{$H+}
{$APPTYPE CONSOLE}
const
l: UInt64 = 6148914691236517205;
var
s,n: UInt64;
i,j: UInt64;
begin
WriteLn('Experiment:', LineEnding);
WriteLn(' Draw a natural number r from the intervall [0,l-1] and');
WriteLn(' increment a counter s when r < l div 2 is satisfied.');
WriteLn(' Repeat this step n times and calculate the ratio s/n.', LineEnding);
WriteLn(' Expected ratio: ', (l div 2)/l:30, LineEnding);
WriteLn('Input size n':16, 'Observed ratio s/n':30);
l := 6148914691236517205;
for j := 4 to 63 do
begin
n := (Int64(1) shl j);
s := 0;
for i := 0 to n-1 do
if Random(Int64(l)) < l div 2 then
Inc(s);
WriteLn( (UInt64(1) shl j):16, s/n:30);
end;
end.
Additional information:
However the 32-bit variant does not display this problem, so one could simply apply the same approach to the 64 bit version. E.g. doing something like this:
{$IFDEF CPUX86_64}
{$ASMMODE INTEL}
function Random64(k: UInt64): UInt64;
var
a: UInt64;
begin
a := (UInt64(GenRand32) shl 32) or GenRand32;
asm
MOV RAX, a
MUL k
MOV result, RDX
end['RAX', 'RDX'];
end;
{$ENDIF}
Or in pure pascal:
function Random64(k: UInt64): UInt64;
var
a, b, c, d: UInt32;
bd, ad, bc, ac: UInt64;
carry: UInt64;
begin
a := GenRand32;
b := GenRand32;
c := k shr 32;
d := UInt32(k);
bd := b * d;
ad := a * d;
bc := b * c;
ac := a * c;
// We only need the carry bit
carry := ((bd shr 32) + UInt32(ad) + UInt32(bc)) shr 32;
// Calculate the final result
Result := carry + (ad shr 32) + (bc shr 32) + ac;
end;
Is there a specific reason the 64 bit variant was implemented using the mod operator?
Mantis conversion info:
- Mantis ID: 35878
- OS: Windows 10 64 bit
- Version: 3.0.4
- Fixed in revision: 42508 (#789f2887)