This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |
Other format: | [Raw text] |
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
Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|
Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |