[Bug-gsl]Bug in discrete random distribution

[ahoward]

Hi folks:

This bug report got bounced from address@hidden, so I'm re-posting it
here...

        Andrew



Hi folks:

I have found a bug in the discrete random number distribution code of your
otherwise excellent package.  For near-uniform distributions,
gsl_ran_discrete() sometimes generates invalid numbers (outside the range
0 <= k < K).  Digging around in the code, I found this:

randist/discrete.c:221

/*** Begin Walker's Algorithm ***/

gsl_ran_discrete_t *
gsl_ran_discrete_preproc(size_t Kevents, const double *ProbArray)

*** snip ***

    /* Now work through the smalls */
    while (size_stack(Smalls) > 0) {
        s = pop_stack(Smalls);
        if (size_stack(Bigs) == 0) {
            /* Then we are on our last value */
            (g->A)[s]=s;
            (g->F)[s]=1.0;
                 break;
        }

*** snip ***

    }
    while (size_stack(Bigs) > 0) {
        b = pop_stack(Bigs);
        (g->A)[b]=b;
        (g->F)[b]=1.0;
    }
    /* Stacks have been emptied, and A and F have been filled */

The problem here is that the Smalls stack may be non-empty on loop
termination, and hence the A and F arrays may be only partially
initialized.  As a work-around, I added:

    while (size_stack(Smalls) > 0) {
        s = pop_stack(Smalls);
        (g->A)[s]=s;
        (g->F)[s]=1.0;
    }

This snippet prevents the generation of garbage, but I cannot speak for
its algorithmic correctness.

Note that I am working with gsl-1.3, but this bug appears in earlier
versions also.

        Cheers,
                Andrew


Andrew Howard                           email: address@hidden
Department of Computer Science          http:  www-robotics.usc.edu/~ahoward
University of Southern California       phone: 1 (213) 740 6416
Los Angeles, CA, U.S.A. 90089-0781      fax:   1 (213) 821 5696
    << Insert pithy saying here >>>