Subtle architecture-dependent bug in Numeric.RandomArray

[Johannes Nix]

There is a subtle bug in the RandomArray module, which is part of
the Numeric package. It causes the random number generator 
RandomArray.normal() to return incorrect values when invoked on Linux 
compiled with gcc  on a AMD Opteron machine, that is, a system with 
64-Bit CPU and ILP64 data model. The result will depend on the data model 
and the aligning of the C implementation on the used architecture. 
The bug seems not to be triggered on 32-bit systems.

Cause of the bug are the code lines

ranlibmodule.c:get_continuous_random, line 48:
---------------------------------------------------
    case 0:
      *out_ptr = (double) ((double (*)()) fun)();
      break;
---------------------------------------------------

and

ranlibmodule.c:standard_normal, line 125:
------------------------------------------------
static PyObject *
standard_normal(PyObject *self, PyObject *args) {
  return get_continuous_random(0, self, args, snorm);
}
-----------------------------------------------------

where fun is a function pointer to the functions ranf() or 
snorm(). However, only in the case of ranf(),
the function returns actually a double value, snorm(),
which is used to generate Gaussian random numbers, returns 
a float.

The Numpy package, which is the recommended replacement
for Numeric, uses a different random number generator
(the Mersenne Twister RNG) and does not contain this code.


Kind Regards,

Johannes Nix