How to implement a fast RNG?

Question

I am trying to port an existing random generator based on 128 bit XorShift from C. But I have trouble with generating the seed which is just generating the same number again and again.

static uint64_t s[ 2 ];

static uint64_t __inline next(void) {
    uint64_t s1 = s[ 0 ];
    const uint64_t s0 = s[ 1 ];
    s[ 0 ] = s0;
    s1 ^= s1 << 23;
    return ( s[ 1 ] = ( s1 ^ s0 ^ ( s1 >> 17 ) ^ ( s0 >> 26 ) ) ) + s0;
}

uint64_t getusertime() {
    struct rusage rusage;
    getrusage( 0, &rusage );
    return rusage.ru_utime.tv_sec * 1000000ULL + ( rusage.ru_utime.tv_usec / 1000 ) * 1000;
}

int main( int argc, char* argv[] ) {
    const long long int n = strtoll( argv[1], NULL, 0 );
    uint64_t t = 0;

    for( int i = 0; i < 2; i++ ) s[ i ] = -1ULL / 3;

    const int64_t start = getusertime();

    for( long long int i = n; i-- != 0; ) t ^= next();

    const int64_t elapsed = getusertime() - start;
    const double secs = elapsed / 1E6;
    printf( "%f s, %.02f queries/s, %.02f ns/query\n", secs, n / secs, 1E9 * secs / n );
    if ( t == 0 ) putchar( 0 );
    return 0;
}

program Project1;

var
  S: Array [0..1] of UInt64;

function XorShift128: UInt64;
var
  s0, s1: UInt64;
begin
  s1 := s[0];
  s0 := s[1];
  s[0] := s0;
  s1 := s1 xor (s1 shl 23);
  s[1] := (s1 xor s0 xor (s1 shr 17) xor (s0 shr 26));
  Result := s[1] + s0;
end;

procedure GenerateSeed;
var
  I: Integer;
begin
  for I := 0 to High(S) do
  S[I] := MaxLongInt div 3;
end;

var
  I: UInt64;
begin 
  GenerateSeed;
  I := XorShift128;
end.

Answer 1

The reason you get the same value every time you run the program in the question is that you use the same seed every time. If I am understanding your comments correctly. The other difference between the C and the Pascal is the seed – see below.

However, your code is fine and is an accurate translation of the C code. The output of this C program:

#include <stdio.h>
#include <stdint.h>

static uint64_t s[ 2 ];

static uint64_t __inline next(void) {
    uint64_t s1 = s[ 0 ];
    const uint64_t s0 = s[ 1 ];
    s[ 0 ] = s0;
    s1 ^= s1 << 23;
    return ( s[ 1 ] = ( s1 ^ s0 ^ ( s1 >> 17 ) ^ ( s0 >> 26 ) ) ) + s0;
}

int main(void) 
{
    s[ 0 ] = s[ 1 ] = 715827882; // the value of MaxLongInt div 3
    printf("%llu\n", next());
    printf("%llu\n", next());
    printf("%llu\n", next());
    return 0;
}

is

6004846026386057
6004846115863870
12676181551404632061

The output of this Delphi program:

program Project1;

{$APPTYPE CONSOLE}

var
  S: Array [0..1] of UInt64;

function XorShift128: UInt64;
var
  s0, s1: UInt64;
begin
  s1 := s[0];
  s0 := s[1];
  s[0] := s0;
  s1 := s1 xor (s1 shl 23);
  s[1] := (s1 xor s0 xor (s1 shr 17) xor (s0 shr 26));
  Result := s[1] + s0;
end;

procedure GenerateSeed;
var
  I: Integer;
begin
  for I := 0 to High(S) do
    S[I] := MaxLongInt div 3;
end;

begin
  GenerateSeed;
  Writeln(XorShift128);
  Writeln(XorShift128);
  Writeln(XorShift128);
end.

is

6004846026386057
6004846115863870
12676181551404632061

I note that the C code in the question uses a different seed from your translation. It seeds the state with -1ULL / 3 and that leads to this output:

46820872945684
46912499612351
13066320939010318272

To match that in the Delphi code you would use high(UInt64) div 3. Do that and you get the output above.

An important note here is that your Delphi code only supplies 64 bits of seed, but your C code supplies 128. I expect that you should supply 128 bits of seed.