bug report on random number generators

bug report on random number generators

[Heiko Bauke]

[-- Attachment #1: Type: text/plain, Size: 4007 bytes --]

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 #2: bug_report_2.tar.gz --]
[-- Type: application/gzip, Size: 3819 bytes --]

searching for gsl library

[Juan Jose Gomez Cadenas]

> Hi,
I have problems to find gsl library using the macro
 AC_CHECK_LIB(library, function)

for instante, I have tried

AC_CHECK_LIB(gsl, gsl_root_fsolver_root)

and many other functions. I have tried in my laptop (gsl installed
in /usr/local/lib ) and at CERN AFS, givin explicitely the gsl path
to LDLFLAGS. It doesn't work in either case.

Instead, with:
AC_CHECK_LIB(gslcblas,cblas_sspmv)

I find gslcblas all right.

Any suggestion?? What am I doing wrong?

I copy configure.in of my application for completness:

dnl Process this file with autoconf to produce a configure script.
AC_PREREQ(2.53)
AC_INIT(version.cpp)
AC_CANONICAL_TARGET([])

dnl Initi automake with the name of the package and the version
AM_INIT_AUTOMAKE(anadiff,"0.1")
AM_CONFIG_HEADER(config.h)

dnl things required by automake

AC_PROG_MAKE_SET([])
AC_PROG_LIBTOOL

dnl Checks for programs.
AC_LANG([C++])
AC_PROG_CC
AC_PROG_CPP
AC_PROG_CXX
AC_PROG_F77

AC_MSG_CHECKING([for math library])
AC_CHECK_LIB(m,main)

AC_MSG_CHECKING([for g2c library])
AC_CHECK_LIB(g2c,main)
AC_PROG_INSTALL
AC_PROG_LN_S
AC_CHECK_TOOL(RANLIB, ranlib, :)
AC_CHECK_TOOL(AR, ar, :)
AC_HEADER_STDC

....  I find OK math and g2c ....

LDFLAGS="$LDFLAGS -L$CERNPATH"
AC_MSG_CHECKING([for packlib library])
AC_CHECK_LIB(packlib,hbook1_)

if test x$GSLPATH = x; then
   GSLPATH=/usr/local/lib
  echo "set GSLPATH to "$GSLPATH
else
     GSLPATH=$GSLPATH
  echo "set GSLPATH to "$GSLPATH
fi

LDFLAGS="$LDFLAGS -L$GSLPATH"

AC_MSG_CHECKING([for gsl library])
AC_CHECK_LIB(gsl, gsl_root_fsolver_root)

I find OK packlib but NOT gsl.

...

 AC_MSG_CHECKING([for gslcblas library])
AC_CHECK_LIB(gslcblas,cblas_sspmv)

LIBS="-lgsl $LIBS"

I find OK gslcblas.

best,

JJ

Re: searching for gsl library

[Martin Jansche]

On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:

> I have problems to find gsl library using the macro

Does the recipe from the manual
(http://sources.redhat.com/gsl/ref/gsl-ref_2.html#SEC12) not work?

- martin

Re: searching for gsl library

[Juan Jose Gomez Cadenas]

> On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:
>
>> I have problems to find gsl library using the macro
>
> Does the recipe from the manual
> (http://sources.redhat.com/gsl/ref/gsl-ref_2.html#SEC12) not work?
>
> - martin

No, it doesn't, at least for me. My configure.in:

AC_MSG_CHECKING([for math library])
AC_CHECK_LIB(m,main)

AC_MSG_CHECKING([for g2c library])
AC_CHECK_LIB(g2c,main)

AC_MSG_CHECKING([for gsl library])
AC_CHECK_LIB(gsl,main)
AM_PATH_GSL
AC_CHECK_LIB(gslcblas,cblas_sdsdot)

ANd the result:
checking for math library... checking for main in -lm... yes
checking for g2c library... checking for main in -lg2c... yes
checking for gsl library... checking for main in -lgsl... no

Re: searching for gsl library

[Alan Aspuru-Guzik]

Hi Juan Jose,

I think you can ommit the AC_CHECK_LIB(gsl,main) completely. AM_PATH_GSL
does that for you...

checking for gsl-config... (cached) /usr/bin/gsl-config
checking for GSL - version >= 0.2.5... yes
checking for inline... (cached) inline

came out from:

dnl Checks for programs.
dnl Checks for libraries.
AM_PATH_GLIB // another lib
AM_PATH_XML2 // another lib
AM_PATH_GSL  // sufficient check for gsl :)

in my configure.in

I also have this line:

AC_CHECK_LIB(atlas,ATL_ctrsv, [ATLAS_LIBS="-latlas
-lcblas"],[ATLAS_LIBS=""])
AC_SUBST(ATLAS_LIBS)

complemented with an $(ATLAS_LIBS) prior to $(GSL_LIBS) in Makefile.am to
automatically check for Atlas.

I think your check fails because there is no symbol named main:

aspuru@cornbread mb]$ nm /usr/lib/libgsl.a | grep main
00002ae0 T gsl_sf_gammainv
00001cc0 T gsl_sf_gammainv_e
         U gsl_sf_gammainv_e
         U gsl_sf_gammainv_e


If you want to use AC_CHECK_LIB, you might use something like:
AC_CHECK_LIB(gsl,gsl_block_alloc) or a function that comes out of 'nm'-ing
the library.

Salúd!

Alán



On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:

> > On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:
> >
> >> I have problems to find gsl library using the macro
> >
> > Does the recipe from the manual
> > (http://sources.redhat.com/gsl/ref/gsl-ref_2.html#SEC12) not work?
> >
> > - martin
> 
> No, it doesn't, at least for me. My configure.in:
> 
> AC_MSG_CHECKING([for math library])
> AC_CHECK_LIB(m,main)
> 
> AC_MSG_CHECKING([for g2c library])
> AC_CHECK_LIB(g2c,main)
> 
> AC_MSG_CHECKING([for gsl library])
> AC_CHECK_LIB(gsl,main)
> AM_PATH_GSL
> AC_CHECK_LIB(gslcblas,cblas_sdsdot)
> 
> ANd the result:
> checking for math library... checking for main in -lm... yes
> checking for g2c library... checking for main in -lg2c... yes
> checking for gsl library... checking for main in -lgsl... no
> 
> 
> 
> 

-- 

Alan Aspuru-Guzik                    Dios mueve al jugador, y éste, la pieza.
(510)642-2154 UC Berkeley           ¿Qué Dios detrás de Dios la trama empieza
(925)422-8739 LLNL                de polvo y tiempo y sueño y agonías? -Borges

Re: searching for gsl library

[Juan Jose Gomez Cadenas]

The following works:

AM_PATH_GSL(0.2.5, LIBS="$LIBS -lgsl")
but one needs the macro gsl.m4 and to specify the version of the library.
I don't understand why I can find gslcblas and not gsl.

Re: searching for gsl library

[Martin Jansche]

On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:

> No, it doesn't, at least for me. My configure.in:
>
> AC_MSG_CHECKING([for math library])
> AC_CHECK_LIB(m,main)
>
> AC_MSG_CHECKING([for gsl library])
> AC_CHECK_LIB(gsl,main)
> AM_PATH_GSL
> AC_CHECK_LIB(gslcblas,cblas_sdsdot)

Try checking for gslcblas before checking for gsl proper.

- martin

Re: searching for gsl library

[Juan Jose Gomez Cadenas]

Martin,
Checking first for gslcblas does it. I imagined that the reason is the
dependency on gsl with gslcblas. Should one perhaps add a little line to
the manual explaining this?

thanks, JJ

> On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:
>
>> No, it doesn't, at least for me. My configure.in:
>>
>> AC_MSG_CHECKING([for math library])
>> AC_CHECK_LIB(m,main)
>>
>> AC_MSG_CHECKING([for gsl library])
>> AC_CHECK_LIB(gsl,main)
>> AM_PATH_GSL
>> AC_CHECK_LIB(gslcblas,cblas_sdsdot)
>
> Try checking for gslcblas before checking for gsl proper.
>
> - martin

Re: searching for gsl library

[Juan Jose Gomez Cadenas]

Thanks Alan, AM_PATH_GSL does the trick. AC_CHECK_LIB works also but one
has to AC_CHECK_LIB firs for gslcblas.

salud,

JJ
> Hi Juan Jose,
>
> I think you can ommit the AC_CHECK_LIB(gsl,main) completely.
> AM_PATH_GSL does that for you...
>
> checking for gsl-config... (cached) /usr/bin/gsl-config
> checking for GSL - version >= 0.2.5... yes
> checking for inline... (cached) inline
>
> came out from:
>
> dnl Checks for programs.
> dnl Checks for libraries.
> AM_PATH_GLIB // another lib
> AM_PATH_XML2 // another lib
> AM_PATH_GSL  // sufficient check for gsl :)
>
> in my configure.in
>
> I also have this line:
>
> AC_CHECK_LIB(atlas,ATL_ctrsv, [ATLAS_LIBS="-latlas
> -lcblas"],[ATLAS_LIBS=""])
> AC_SUBST(ATLAS_LIBS)
>
> complemented with an $(ATLAS_LIBS) prior to $(GSL_LIBS) in Makefile.am
> to automatically check for Atlas.
>
> I think your check fails because there is no symbol named main:
>
> aspuru@cornbread mb]$ nm /usr/lib/libgsl.a | grep main
> 00002ae0 T gsl_sf_gammainv
> 00001cc0 T gsl_sf_gammainv_e
>         U gsl_sf_gammainv_e
>         U gsl_sf_gammainv_e
>
>
> If you want to use AC_CHECK_LIB, you might use something like:
> AC_CHECK_LIB(gsl,gsl_block_alloc) or a function that comes out of
> 'nm'-ing the library.
>
> Salúd!
>
> Alán
>
>
>
> On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:
>
>> > On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:
>> >
>> >> I have problems to find gsl library using the macro
>> >
>> > Does the recipe from the manual
>> > (http://sources.redhat.com/gsl/ref/gsl-ref_2.html#SEC12) not work?
>> >
>> > - martin
>>
>> No, it doesn't, at least for me. My configure.in:
>>
>> AC_MSG_CHECKING([for math library])
>> AC_CHECK_LIB(m,main)
>>
>> AC_MSG_CHECKING([for g2c library])
>> AC_CHECK_LIB(g2c,main)
>>
>> AC_MSG_CHECKING([for gsl library])
>> AC_CHECK_LIB(gsl,main)
>> AM_PATH_GSL
>> AC_CHECK_LIB(gslcblas,cblas_sdsdot)
>>
>> ANd the result:
>> checking for math library... checking for main in -lm... yes
>> checking for g2c library... checking for main in -lg2c... yes
>> checking for gsl library... checking for main in -lgsl... no
>>
>>
>>
>>
>
> --
>
> Alan Aspuru-Guzik                    Dios mueve al jugador, y éste, la
> pieza. (510)642-2154 UC Berkeley           ¿Qué Dios detrás de Dios la
> trama empieza (925)422-8739 LLNL                de polvo y tiempo y
> sueño y agonías? -Borges

Re: searching for gsl library

[Martin Jansche]

On Fri, 21 Feb 2003, Juan Jose Gomez Cadenas wrote:

> Checking first for gslcblas does it. I imagined that the reason is the
> dependency on gsl with gslcblas. Should one perhaps add a little line to
> the manual explaining this?

TBH, I think it may already be clear from the manual.  But then again,
you've just demonstrated that some further clarification is desirable.

- martin

Re: bug report on random number generators

[Brian Gough]

Heiko Bauke writes:
 > 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.
 > 

Thanks for the bug report and patches. 

Sorry for not replying to your previous message.  I forwarded it to
the original author of the Knuth routines, Carlo Perassi, and he is
working on fixing them.  

Ideally these bugs should have been caught before they were included
by the test of computing the 10,000th value analytically but I must
have let them slip through without checking that.

I will correct the period in the comments from 2^32 to 2^30.  Thanks
for pointing that out.  

I don't use (s|1) since it produces identical sequences for even and
odd seeds.  Although this allows seeds with sub-maximal period this is
the intention -- the 'randu' type generators are provided as a way to
produce poor-quality sequences.

-- 
Brian

bug in covariance_source.cpp??

[Juan Jose Gomez Cadenas]

In covariance_source.cpp, the function:

double
FUNCTION(gsl_stats,covariance) (const BASE data1[], const size_t stride1,
                                const BASE data2[], const size_t stride2,
                                const size_t n)
{
  const double mean1 = FUNCTION(gsl_stats,mean) (data1, stride1, n);
  const double mean2 = FUNCTION(gsl_stats,mean) (data2, stride2, n);

  return FUNCTION(gsl_stats,covariance_m)(data1, stride1,
                                          data2, stride2,
                                          n,
                                          mean1, mean2);
}

calls the function covariance_m:

double
FUNCTION(gsl_stats,covariance_m) (const BASE data1[], const size_t
stride1,
                                  const BASE data2[], const size_t
stride2,
                                  const size_t n,
                                  const double mean1, const double mean2)
{
  const double covariance = FUNCTION(compute,covariance) (data1, stride1,
                                                          data2, stride2,
                                                          n,
                                                          mean1, mean2);

  return covariance * ((double)n / (double)(n - 1));
}


but I believe it should call the function covariance:
static double
FUNCTION(compute,covariance) (const BASE data1[], const size_t stride1,
                              const BASE data2[], const size_t stride2,
                              const size_t n,
                              const double mean1, const double mean2)
{
  /* takes a dataset and finds the covariance */

  long double covariance = 0 ;

  size_t i;

  /* find the sum of the squares */
  for (i = 0; i < n; i++)
    {
      const long double delta1 = (data1[i * stride1] - mean1);
      const long double delta2 = (data2[i * stride2] - mean2);
      covariance += (delta1 * delta2 - covariance) / (i + 1);
    }

  return covariance ;
}

Regards, JJ

Re: bug in covariance_source.cpp??

[Brian Gough]

Juan Jose Gomez Cadenas writes:
 > In covariance_source.cpp, the function gsl_stats_covariance
 > calls the function covariance_m but I believe it should
 > call the function covariance


The _m function takes a precomputed mean, rather than a "known
mean" (in the statistical sense), so I believe the factor of
1/(n-1) is ok and should be the same for both cases.  Maybe
this should be clarified in the documentation though.

-- 
Brian

Re: bug report on random number generators

[Brian Gough]

Heiko Bauke writes:
 > 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.

These are now fixed in CVS ready for the next release.

-- 
Brian