[Bug-gsl] a seeming bug with gsl_ran_gaussian and similar functions

[Vladimir L. Koliadin]

=================
1. gsl-1.8.tar.gz  from www.gnu.org 
      (and the same bug has been observed 
       with the gsl-1.7 included into FC5 distribution)

=================
2. Hardware:  Pentium 4 3GHz, Intel D945 PSNLK, 1GB RAM 
   (see details below)

   Operating system: Fedora Core 5 (FC5) 
   
uname -a
>   Linux pc2 2.6.17-vk #1 SMP Sun Aug 13 21:42:26 EEST 2006 
>    i686 i686 i386 GNU/Linux


================
3. Compiler: 

gcc --version
> gcc (GCC) 4.1.0 20060304 (Red Hat 4.1.0-3)

Compiler options (two variants):
   CFLAGS=" -g -O2"
   CFLAGS=" -g "

================
4. Bug behavior: 

Functions ``gsl_ran_gaussian'' and ``gsl_ran_exponential'' 
produce definitely wrong results if called directly (i.e. from the GSL library)

But, if the code of the function is copied (under another function name) 
to the test program, then the results are OK.

The bug seems to be related to gsl_rng_uniform_pos():
at least it is not observed with functions like gsl_ran_poisson(),
which don't call gsl_rng_uniform_pos().


===============
5. Bug-demo program:

---------------- begin prgram ---------------------------------
/*
   FILE: gsl-bug.c

   demostration of a (seeming?) bug with functions like:
      gsl_ran_gaussian
      gsl_ran_exponential
      ...

   the bug seems to be related to gsl_rng_uniform_pos(...)
   (because it is not observed with functions like gsl_ran_poisson(...),
    which don't refer to gsl_rng_uniform_pos(...) )

   COMPILE: cc -o a -g gsl-bug.c -lgsl -lgslcblas -lm 

   try it with ``use_copy = 1'' (all is OK) and ``use_copy=0''

---   The output with use_copy = 1: (OK) ---
0  1.085807
1  -0.566087
2  0.729758
3  1.739741
4  -0.073342

---   The output with use_copy = 0: (WRONG) ---
0  0.000000
1  0.000000
2  0.000000
3  0.000000
4  0.000000

*/

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <math.h>

int use_copy = 0;  // if 1, then the gsl_ran_gaussian_copy defined below 
                   // if 0, then the gsl_ran_gaussian itself is used

// it is a (renamed) copy of the function ``gsl_ran_gaussian(...)''  from
//    the source distribution of gsl-1.8
double gsl_ran_gaussian_copy (const gsl_rng * r, const double sigma)
{
  double x, y, r2;

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
      x = -1 + 2 * gsl_rng_uniform_pos (r);
      y = -1 + 2 * gsl_rng_uniform_pos (r);

      /* see if it is in the unit circle */
      r2 = x * x + y * y;
    }
  while (r2 > 1.0 || r2 == 0);

  /* Box-Muller transform */
  return sigma * y * sqrt (-2.0 * log (r2) / r2);
}

int main (void)
{
  const gsl_rng_type * T = gsl_rng_mt19937; // allocate (create)  rng:
  gsl_rng * r = gsl_rng_alloc(T);           // specify rng-type:
  unsigned long int seed=123;
  int i, n = 5; 

  gsl_rng_set(r, seed);  // set seed 

  if(use_copy) // call the copy of the function
  {
    for(i=0; i<n; i++) 
         printf("%d  %f\n", i, gsl_ran_gaussian_copy(r, 1.0));
  }
  else // call the original function from the GLS library
  {
    for(i=0; i<n; i++) 
         printf("%d  %f\n", i, gsl_ran_gaussian(r, 1.0));
  }
};

---------------- end prgram ---------------------------------



=======================
Hardware details:

======
cat /proc/cpuinfo

> processor     : 0
> vendor_id     : GenuineIntel
> cpu family    : 15
> model         : 4
> model name    : Intel(R) Pentium(R) 4 CPU 3.00GHz
> stepping      : 3
> cpu MHz               : 3000.821
> cache size    : 2048 KB
> physical id   : 0
> siblings      : 2
> core id               : 0
> cpu cores     : 1
> fdiv_bug      : no
> hlt_bug               : no
> f00f_bug      : no
> coma_bug      : no
> fpu           : yes
> fpu_exception : yes
> cpuid level   : 5
> wp            : yes
> flags         : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov 
> pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe nx lm constant_tsc 
> pni monitor ds_cpl est cid cx16 xtpr
> bogomips      : 6008.26
> 
> processor     : 1
> vendor_id     : GenuineIntel
> cpu family    : 15
> model         : 4
> model name    : Intel(R) Pentium(R) 4 CPU 3.00GHz
> stepping      : 3
> cpu MHz               : 3000.821
> cache size    : 2048 KB
> physical id   : 0
> siblings      : 2
> core id               : 0
> cpu cores     : 1
> fdiv_bug      : no
> hlt_bug               : no
> f00f_bug      : no
> coma_bug      : no
> fpu           : yes
> fpu_exception : yes
> cpuid level   : 5
> wp            : yes
> flags         : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov 
> pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe nx lm constant_tsc 
> pni monitor ds_cpl est cid cx16 xtpr
> bogomips      : 6000.58


======
/sbin/lspci

> 00:00.0 Host bridge: Intel Corporation 945G/P Memory Controller Hub
> 00:01.0 PCI bridge: Intel Corporation 945G/P PCI Express Graphics Port
> 00:1b.0 Audio device: Intel Corporation 82801G (ICH7 Family) High Definition 
> Audio Controller (rev 01)
> 00:1c.0 PCI bridge: Intel Corporation 82801G (ICH7 Family) PCI Express Port 1 
> (rev 01)
> 00:1c.2 PCI bridge: Intel Corporation 82801G (ICH7 Family) PCI Express Port 3 
> (rev 01)
> 00:1c.3 PCI bridge: Intel Corporation 82801G (ICH7 Family) PCI Express Port 4 
> (rev 01)
> 00:1d.0 USB Controller: Intel Corporation 82801G (ICH7 Family) USB UHCI #1 
> (rev 01)
> 00:1d.1 USB Controller: Intel Corporation 82801G (ICH7 Family) USB UHCI #2 
> (rev 01)
> 00:1d.2 USB Controller: Intel Corporation 82801G (ICH7 Family) USB UHCI #3 
> (rev 01)
> 00:1d.3 USB Controller: Intel Corporation 82801G (ICH7 Family) USB UHCI #4 
> (rev 01)
> 00:1d.7 USB Controller: Intel Corporation 82801G (ICH7 Family) USB2 EHCI 
> Controller (rev 01)
> 00:1e.0 PCI bridge: Intel Corporation 82801 PCI Bridge (rev e1)
> 00:1f.0 ISA bridge: Intel Corporation 82801GB/GR (ICH7 Family) LPC Interface 
> Bridge (rev 01)
> 00:1f.1 IDE interface: Intel Corporation 82801G (ICH7 Family) IDE Controller 
> (rev 01)
> 00:1f.2 IDE interface: Intel Corporation 82801GB/GR/GH (ICH7 Family) Serial 
> ATA Storage Controllers cc=IDE (rev 01)
> 00:1f.3 SMBus: Intel Corporation 82801G (ICH7 Family) SMBus Controller (rev 
> 01)
> 01:00.0 VGA compatible controller: nVidia Corporation NV43 [GeForce 6600] 
> (rev a2)
> 02:00.0 Ethernet controller: Intel Corporation 82573V Gigabit Ethernet 
> Controller (Copper) (rev 03)
> 05:03.0 Multimedia controller: Philips Semiconductors SAA7133 Video Broadcast 
> Decoder (rev f0)
> 05:05.0 FireWire (IEEE 1394): Texas Instruments TSB43AB23 IEEE-1394a-2000 
> Controller (PHY/Link)