bug report on random number generators

[Heiko Bauke <Heiko dot Bauke at Physik dot Uni-Magdeburg dot DE>]

Hi,

I have found some bugs in the initialization routine for generators
gsl_rng_transputer, gsl_rng_randu, gsl_rng_borosh13 and some other
serious bugs.

In general generators of type

x_i = a*x_{i-1} mod 2^e

have a maximal period of 2^{e-2}. A generator has this period if and
only if

* e>2
* +/- 3 = a mod 8  and 
* x_1 is odd.

For generator gsl_rng_transputer e=32, so the period is only 2^30, not
2^32 as the comments in the sources gsl-1.3/rng/transputer.c suggest.

Let's consider generator

      x_i = 11 x_{i-1} mod 2^5 .

With x_0=1 we get the full period 1, 11, 25, 19, 17, 27, 9, 3, 1, ...
but with x_0=12 the sequence jumps between 12 and 4 and we get the
sequence 12, 4, 12, 4, 12,... Note that in a maximal period sequence

* the lowest bit is constant,
* the 2nd lowest bit has a period of 2^1,
* the 3rd lowest is constant,
* the 4th lowest bit has a period of 2^2,
* the 5th lowest bit has a period of 2^3 and so on.

To ensure that the generator has maximal period the initialization
routine 

static void
transputer_set (void *vstate, unsigned long int s)
{
  transputer_state_t *state = (transputer_state_t *) vstate;

  if (s == 0)
    s = 1 ;   /* default seed is 1. */

  state->x = s;

  return;
}


should be changed into


static void
transputer_set (void *vstate, unsigned long int s)
{
  transputer_state_t *state = (transputer_state_t *) vstate;

  state->x = s|1;  /* odd seed */

  return;
}


I also found bugs in the linear congruential generators with modulus 
2^31-1. These generators actually calculate modulo 2^31. That's
something completely different. This bug concerns to fishman18,
fishman20, fishman2x and kunthran2.

The implementations of the initialization routines of some linear
congruential generators are also problematic. 

I would like to contribute some bug fixes. Here a list of changes I
have made:

* borosh13
  - initialization routine changed to ensure a odd seed
  - RAND_MIN = 1

* coveyou
  - RAND_MAX = 2^32-2
  - RAND_MIN = 2
  ( coveyou generates only numbers x_i with  2 = x_i mod 4 )

* fishman18
 - RAND_MAX = 2^31-2
 - RAND_MIN = 1
 - calculates now   x_i = 62089911 x_{i-1} mod (2^31-1)
   ( not x_i = 62089911 x_{i-1} mod (2^31) ) as described by Knuth
 - initialization routine changed to ensure 0 < x < (2^31-1)
 - ran_get_double modified

* fishman20
 - RAND_MAX = 2^31-2
 - RAND_MIN = 1
 - calculates now   x_i = 48271 x_{i-1} mod (2^31-1)
   ( not x_i = 48271 x_{i-1} mod (2^31) ) as described by Knuth
 - initialization routine changed to ensure 0 < x < (2^31-1)
 - ran_get_double modified

* fishman2x
 - RAND_MAX = 2^31-2
 - calculates now   
        x_i = 48271 x_{i-1} mod (2^31-1)
        y_i = 40692 y_{i-1} mod (2^31-249)
        z_i = x_i - y_i mod (2^31-1)
   (not  x_i = 48271 x_{i-1} mod (2^31), y_i = 40692 y_{i-1} mod (2^31-249),
     x_i = x_i - y_i mod (2^31) ) as described by Knuth
 - initialization routine changed to ensure 0 < x < (2^31-1) and 
   0 < y < (2^31-249) 
 - ran_get_double modified 

* kunthran2
 - RAND_MAX = 2^31-2
 - calculates now  x_i = 271828183 x_{i-1} - 314159269 x_{i-2} mod (2^31-1)
   ( not x_i = 271828183 x_{i-1} - 314159269 x_{i-2} mod (2^31) ) as
   described by Knuth 
 - initialization routine changed to ensure 0 < x0 < (2^31-1) and 
   0 < x1 < (2^31-1) 
 - ran_get_double modified 

* minst
 - initialization routine changed to ensure 0 < x < (2^31-1)

* lecuyer21
 - RAND_MAX = 2^31-250
 - RAND_MIN = 1
 - initialization routine changed to ensure 0 < x < (2^31-249)
 - ran_get_double modified 

* transputer
 - initialization routine changed to ensure a odd seed

* waterman
 - RAND_MIN = 1
 - initialization routine changed to ensure a odd seed


Remark on mrg:
The quality of mrg does not depend on the 5 inital values as long not
all x_1, ... , x_5 are zero but smaller than 2^31-1. The warm up
procedure is not needed. 


	Heiko



-- 
-- Supercomputing in Magdeburg @ http://tina.nat.uni-magdeburg.de
--                 Heiko Bauke @ http://www.uni-magdeburg.de/bauke

Attachment: bug_report_2.tar.gz
Description: application/gzip