[Help-gsl] Possible bug with the "gsl_rng_ranlxs2" routine

[Wang, Gao]

Hi,

I am using GSL-1.14 and I got strange result when I generate pseudo random
numbers X~U(0, 1). The problem occurs randomly at different runs of the
problem. In example A below it is problematic, and B looks fine. My test C++
code w/ a make file is attached at the end of the email.

Any comments? Thank you!

### Example A####
./ranlxsbug
use the gsl_rng_ranlxs2 method
Myseed = 3097929819
myunif = 0
myunif = -5.004904628e+12
myunif = 0
myunif = 5.886644847e+12
myunif = 0
myunif = 9.517631799e+12
myunif = 0
myunif = -1.157594218e+12
myunif = 0
myunif = -1.197020912e+13
myunif = 0
myunif = -5.934653374e+12
myunif = 0
myunif = 1.05405093e+13
myunif = 0
use the gsl_rng_mt19937 method
Myseed = 3097929819
myunif = 0.3331135246
myunif = 0.4753004052
myunif = 0.6840925638
myunif = 0.6905622678
myunif = 0.4087153792
myunif = 0.4675788379
myunif = 0.06461505266
myunif = 0.06030185358
myunif = 0.5714081186
myunif = 0.7716544913
myunif = 0.1229928406
myunif = 0.09579540906
myunif = 0.7758432983
myunif = 0.432565941
myunif = 0.8234983671
### Example B####
./ranlxsbug
use the gsl_rng_ranlxs2 method
Myseed = 106374854
myunif = 0.6294272542
myunif = 0.02813184261
myunif = 0.3676368594
myunif = 0.2936396003
myunif = 0.857829988
myunif = 0.9503712058
myunif = 0.6056731343
myunif = 0.1007401943
myunif = 0.4769713283
myunif = 0.3925215602
myunif = 0.7966954112
myunif = 0.9465102553
myunif = 0.4553982019
myunif = 0.3809945583
myunif = 0.6535906792
use the gsl_rng_mt19937 method
Myseed = 106374854
myunif = 0.3676483827
myunif = 0.04865595303
myunif = 0.7878059589
myunif = 0.7906315061
myunif = 0.1233621615
myunif = 0.1731820586
myunif = 0.8502643262
myunif = 0.4718055225
myunif = 0.4755936258
myunif = 0.9699133912
myunif = 0.5249320145
myunif = 0.9644626162
myunif = 0.9786006608
myunif = 0.4194381249
myunif = 0.7208476516
#######################

//!\file BUG_gsl_rng_ranlxs2.cpp
#include <iostream>
#include <cstdlib>
#include <cstddef>
#include <ctime>
#include "gsl/gsl_rng.h"
gsl_rng* initialize_gslr(int specifyRng)
{
  gsl_rng* gslr;
  if (specifyRng == 0)
      gslr = gsl_rng_alloc(gsl_rng_ranlxs2);
  else
      gslr = gsl_rng_alloc(gsl_rng_mt19937);
  unsigned long seed = static_cast<unsigned long>(time (NULL) * getpid() +
10);
  std::cout << "Myseed = " << seed << std::endl;
  gsl_rng_set (gslr, seed);
  return gslr;
}

int main()
{
  std::cout.precision(10);

  std::cout << "use the gsl_rng_ranlxs2 method" << std::endl;
  gsl_rng* gslr = initialize_gslr(0);

  for (int i = 0; i != 15; ++i) {
    double runif = gsl_rng_uniform(gslr);
    std::cout << "myunif = " << runif << std::endl;
  }
  gsl_rng_free(gslr);

  std::cout << "use the gsl_rng_mt19937 method" << std::endl;
  gslr = initialize_gslr(1);

  for (int i = 0; i != 15; ++i) {
    double runif = gsl_rng_uniform(gslr);
    std::cout << "myunif = " << runif << std::endl;
  }
  gsl_rng_free(gslr);

  return 0;
}

//!\file makefile
EXE:= ranlxsbug

SRC:= BUG_gsl_rng_ranlxs2.cpp
OBJS:= $(SRC:.cpp=.o)

CC:= gcc        # The C compiler.
CFLAGS:= -O2 -Wall    # C compilation options which relate to
                # optimization or debugging (usually
                # just -g or -O).  Usually this wouldn't
                # include -I options to specify the
                # include directories, because then you
                # couldn't override it on the command line
                # easily as in the above example.
CXX:= g++         # The C++ compiler.  (Sometimes "CPP" instead
                # of CXX.)
CXXFLAGS:= -O2 -Wall    # C++ compilation options related to
                # optimization or debugging (-O or -g).



INCLUDEDIR:= -I/usr/local/include
LIBDIR:= -L/usr/local/lib
                # Contains libraries we need to link in.
LIBS:= -lgsl -lgslcblas -lm

help :
    @echo " "
    @echo "Type ...          To ..."
    @echo "make all          Compile the program"
    @echo "make clean        Delete temporary files"
    @echo " "

all: $(EXE)

$(EXE): $(OBJS)
    $(CXX) $(INCLUDEDIR) $(OBJS) -o $(EXE) $(LIBDIR) $(LIBS)

.cpp.o: $*.cpp
    $(CXX) $(CXXFLAGS) $(INCLUDEDIR) -o $@ -c $*.cpp
.SUFFIXES : .cpp .c .o $(SUFFIXES)

clean:
    rm *.o
############################


Kindest regards,
Gao

Student in Statistical Genetics, Baylor College of Medicine
(the same Gao as from: gaow [at] bcm.edu / gaow [at] rice.edu)