[GRASS-dev] r.mapcalc random result different on MS Windows with same seed

[Glynn Clements]

Vaclav Petras wrote:

> it seems that r.mapcalc gives different (pseudo-)random result on MS
> Windows then it is giving on Linux. In NC sample dataset:
> 
> g.region n=20 s=10 e=25 w=15 res=1
> r.mapcalc "rand_cell = rand(1, 200)" seed=500
> r.info rand_cell -r
> 
> gives the following on Linux:
> 
> min=1
> max=195
> 
> while on MS Windows it gives:
> 
> min=1
> max=196

I believe that this is caused by (64-bit) Linux having a 64-bit "long"
while Windows (32-bit and 64-bit) has a 32-bit "long".

r.mapcalc's rand() function (f_rand(), in xrand.c) does:

	unsigned long x = (unsigned long)G_mrand48();
	int lo = arg1[i];
	int hi = arg2[i];

	if (lo > hi) {
	    int tmp = lo;

	    lo = hi;
	    hi = tmp;
	}
	res[i] = (lo == hi) ? lo : lo + x % (unsigned long)(hi - lo);

G_mrand48() returns a signed 32-bit result in the range [-2^31, 2^31). 
If the result is negative, casting to an unsigned 32-bit type will
yield a value in the range 0x80000000..0xFFFFFFFF while casting to an
unsigned 64-bit type will yield a value in the range
0xFFFFFFFF80000000..0xFFFFFFFFFFFFFFFF. Unless hi-lo is a power of
two, the result modulo (hi-lo) will be different depending upon
whether the intermediate value is 32-bit or 64-bit.

r64709 replaces the use of "unsigned long" with "unsigned int", which
will be 32-bit on all common platforms. This results in the above test
producing max=196 on (64-bit) Linux.

Technically, this still isn't portable, as converting a negative
integer to an unsigned type invokes undefined behaviour according to
the C standards. In practice, the results will be consistent on any
architecture which has a 32-bit "int" and uses two's complement
representation for signed values, which is probably everything on
which GRASS will run (and anything else is likely to encounter far
bigger issues than r.mapcalc's rand() function).

An alternative fix would be to use G_lrand48() instead of G_mrand48(). 
The former returns a non-negative result in the range [0, 2^31), i.e. 
there are only 2^31 possible values rather than 2^32.

Note that G_lrand48() and G_mrand48() have "long" as their return type
to match the POSIX lrand48() and mrand48() functions (with which they
should be compatible; for the same seed, the GRASS versions should
return the same result as the POSIX versions). Early Unix systems
often had a 16-bit "int", hence the use of "long" in spite of the fact
that the results are constrained to a signed 32-bit range.

-- 
Glynn Clements <glynn at gclements.plus.com>