random numbers

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

Hi,

I have found some bugs in GSL related to pseudo random number generation.

In the GSL documentation one can read: (generator gsl_rng_gfsr4)

"If the offsets are appropriately chosen (such the one ones in this
implementation), then the sequence is said to be maximal. I'm not sure
what that means, but I would guess that means all states are part of
the same cycle, which would mean that the period for this generator is
astronomical; it is @c{$(2^K)^D \approx 10^{93334}$} @math{(2^K)^D
\approx 10^@{93334@}} where @math{K=32} is the number of bits in the
word, and D is the longest lag.  This would also mean that any one
random number could  easily be zero; ie @c{$0 \le r < 2^{32}$} 
@math{0 <= r < 2^32}."

That's not right. The period is much shorter.

"Ziff doesn't say so, but it seems to me that the bits are completely
independent here, so one could use this as an efficient bit generator;
each number supplying 32 random bits."

Yes, this generator is indeed just a combination of 32 one bit
generators. For this reason the period is "only" 2^9689 - 1, but it
has a loot of independent cycles. The cycle is choosen at
initialization time. The period is 2^9689 - 1 and not 2^9689 because
if a seed table with all bits 0 is choosen, the generator produces
only 0s. 

A second error can be found in the documentation for generator
gsl_rng_knuthran2. Parameter a_2 is (according to Knuth) 
a_2 = -314159269 (not 314159269). However in file rng/knuthran2.c a
completely different generator is implementet. Kuth's generator is

  x_n = (271828183*x_{n-1} - 314159269*x_{n-2}) mod m

with m=2^31 - 1. It has a period of m^2 - 1. The generator implemented
in rng/knuthran2.c uses m=2^31. This is a completely different
generator. For example its period is 2*m (or less), see Kuth section
3.2.2 exercise 12. The reason for the different periods is that 2^31-1
is a prime while 2^31 not.

I'd like to contribute a correct (but optimized) implementation of
knuthran2. 

    regards

    Heiko


--- source code for rng/knuthran2.c starts here ---
/* rng/knuthran2.c
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 * This generator is taken from
 *
 * Donald E. Knuth
 * The Art of Computer Programming
 * Volume 2
 * Third Edition
 * Addison-Wesley
 * Page 108
 *
 * This implementation  copyright (C) 2001 Carlo Perassi
 * and (C) 2003 Heiko Bauke.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#define AA1      271828183UL
#define AA2     1833324378UL      /* = -314159269 mod (2 ^ 31 -1) */
#define MM      0x7fffffffUL	  /* 2 ^ 31 - 1 */
#define CEIL_SQRT_MM 46341UL      /* sqrt(2 ^ 31 - 1) */

static inline unsigned long int gsl_knuthran2_schrage(unsigned long int,
unsigned long int);
static inline unsigned long int gsl_knuthran2_mult(unsigned long int, unsigned
long int);
static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);

typedef struct
{
  unsigned long int x0;
  unsigned long int x1;
}
ran_state_t;

static inline unsigned long int
gsl_knuthran2_schrage (unsigned long int a, unsigned long int b) 
{
  /* This is a modified version of Schrage's method. It ensures that no
   * overflow or underflow occurs even if a=ceil(sqrt(MM)). Usual 
   * Schrage's method works only until a=floor(sqrt(MM)).
   */
  unsigned long int q, t;
  if ( a==0UL )
    return 0UL;
  q = MM/a;
  t = 2*MM-(MM%a)*(b/q);
  if (t>=MM)
    t -= MM;
  t += a*(b%q);
  return (t>=MM) ? (t-MM) : t;
}

static inline unsigned long int
gsl_knuthran2_mult (unsigned long int a, unsigned long int b) 
{
  /* To multiply a and b use Schrage's method 3 times.
   * represent a in base ceil(sqrt(m))  a = a1*ceil(sqrt(m)) + a0  
   * a*b = (a1*ceil(sqrt(m)) + a0)*b = a1*ceil(sqrt(m))*b + a0*b   
   */
  unsigned long int t =
    gsl_knuthran2_schrage(a/CEIL_SQRT_MM, 
			  gsl_knuthran2_schrage(CEIL_SQRT_MM, b)) +
    gsl_knuthran2_schrage(a%CEIL_SQRT_MM, b);
  return (t>=MM) ? (t-MM) : t;
}

static inline unsigned long int
ran_get (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;
  
  unsigned long int xtmp = state->x1;
  state->x1 = gsl_knuthran2_mult(AA1, state->x1) + 
    gsl_knuthran2_mult(AA2, state->x0);
  if ( state->x1>=MM )
    state->x1 -= MM;
  state->x0 = xtmp;
  
  return state->x1;
}

static double
ran_get_double (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 2147483647.0;
}

static void
ran_set (void *vstate, unsigned long int s)
{
  ran_state_t *state = (ran_state_t *) vstate;

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

  state->x0 = s % MM;
  state->x1 = s % MM;

  return;
}

static const gsl_rng_type ran_type = {
  "knuthran2",			/* name */
  MM-1L,       			/* RAND_MAX */
  0,				/* RAND_MIN */
  sizeof (ran_state_t),
  &ran_set,
  &ran_get,
  &ran_get_double
};

const gsl_rng_type *gsl_rng_knuthran2 = &ran_type;

--- source code for rng/knuthran2.c ends here ---


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