BUG Report for QRNG-Niederreiter-sequence

BUG Report for QRNG-Niederreiter-sequence

[Wolfgang Hoermann]

/*
This is a bug report and bug-demonstration file for gsl-1.2

Hardware: Old PC
SYSTEM RedHat Linux 6.2

Compiler: gcc (sorry I was not able to determine the version number.
               I installed it directly from the red-hat 6.2 CD) 

compiling command used: gcc -static -Wall testbug.c -lgsl -lm 



Bug description: 
For qrng for the Niedereiter-sequence.
The below program should generate the same sequence 3 times
and thus also print three identical lines, as the state of
the sequence is reset by the   "gsl_qrng_init(qmc)" - command
between the for-loops.

But the output on my computer is

0.171875,0.890625,0.628906,0.378906,
0.171875,0.890625,0.628906,0.378906,
0.171875,0.890625,0.519531,0.378906,

There occur similar problems for dimension 3 and for higher dimensions.

******************************************/

#include <stdio.h>
#include <gsl/gsl_qrng.h>
#define M 40
#define DIM 4

int main (void)
{
  double v[DIM];
  int i;
  gsl_qrng *qmc;

  qmc = gsl_qrng_alloc (gsl_qrng_niederreiter_2, DIM);
  
  for(i=0;i<M;i++) gsl_qrng_get(qmc, v);
  for(i=0;i<DIM;i++) printf("%f,",v[i]);
  printf("\n");

  gsl_qrng_init(qmc);/*resets the state of the "generator" to the
start*/

  for(i=0;i<M;i++) gsl_qrng_get(qmc, v);
  for(i=0;i<DIM;i++) printf("%f,",v[i]);
  printf("\n");

  gsl_qrng_init(qmc);/*resets the state of the "generator" to the
start*/

  for(i=0;i<M;i++) gsl_qrng_get(qmc, v);
  for(i=0;i<DIM;i++)  printf("%f,",v[i]);
  printf("\n");

  gsl_qrng_free(qmc);
  return 0;
}

Re: BUG Report for QRNG-Niederreiter-sequence

[Brian Gough]

Wolfgang Hoermann writes:
 > This is a bug report and bug-demonstration file for gsl-1.2 For
 > qrng for the Niedereiter-sequence.  The below program should
 > generate the same sequence 3 times and thus also print three
 > identical lines, as the state of the sequence is reset by the
 > "gsl_qrng_init(qmc)" - command between the for-loops.
 >  But the output on my computer is
 >  0.171875,0.890625,0.628906,0.378906,
 > 0.171875,0.890625,0.628906,0.378906,
 > 0.171875,0.890625,0.519531,0.378906,
 >  There occur similar problems for dimension 3 and for higher
 > dimensions.

Thanks for the bug report. I found an uninitialised memory access with
checkergcc which looks like to source of the problem.  Here is a
patch.

Index: niederreiter-2.c
===================================================================
RCS file: /cvs/gsl/gsl/qrng/niederreiter-2.c,v
retrieving revision 1.11
diff -u -5 -r1.11 niederreiter-2.c
--- niederreiter-2.c    19 Nov 2001 22:32:08 -0000      1.11
+++ niederreiter-2.c    16 Nov 2002 19:00:57 -0000
@@ -249,10 +249,16 @@
 
     for(k=0; k<=px_degree; k++) {
       px[k] = primitive_poly[poly_index][k];
       pb[k] = 0;
     }
+
+    for (;k<NIED2_MAX_DEGREE+1;k++) {
+      px[k] = 0;
+      pb[k] = 0;
+    }
+
     pb[0] = 1;
 
     for(j=0; j<NIED2_NBITS; j++) {
 
       /* If U = 0, we need to set B to the next power of PX