new randist/binomial.c

[James Theiler <jt at lanl dot gov>]

Hello,
  Awhile back I wrote about how much slower the GSL binomial random
variates were than competing algorithms -- thanks again to John
Pearson for first pointing this out to me.
  I have attached new code for binomial variates which incorporates a
state of the art algorithm (called BTPE [*], the same basic algorithm
that is implemented in ranlib, and I am grateful to Dr.
Kachitvichyanukul and Dr. Schmeiser for permission to modify that code
and include it in GSL; also for useful comments and discussion about
the code and random variate generation in general).  The new code is
about an order of magnitude faster than our original algorithm, and
several times faster than the Numerical Recipes algorithm.  Although
it is essentially the same algorithm that appears in ranlib, I've
managed through various tweaks to produce something that is a little
bit faster than the ranlib code (at least on my machine).
  More details are available in the introductory comments to the
source code, but a few points of general interest:

  1. Following GSL conventions, this algorithm does not employ static
variables.  However, by using static variables to hold intermediate
computations, you can produce significantly (~50%) faster variates in
the situation where you call the function multiple times with the same
parameters.  I have to admit that I found this performance gain too
seductive to resist, so I have included in the source code a '#define
USE_XSTATIC' (the XSTATIC isn't supposed to convey my excitement, just
wanted to avoid possible name conflict with USE_STATIC) which is by
default set to zero, but which an individual can set to 1 to get
faster performance.

  2. As if that wasn't enough, I also implemented another approach --
one that doesn't use static variables -- for speeding up multiple
calls with the same parameters (for binomial, they are n and p).  The
formulation follows that of gsl_ran_discrete (and in fact uses the
gsl_ran_discrete algorithm as well);  take a look at the file
binomial_batch.c for details.

  3. I have added a new function gsl_ran_binomial_cdf(), which is the
cumulative distribution function (akin to the pdf).  I think we should
begin filling out the cdf's for all our nonuniform variates.

  4. Just a quick note, I got a substantial performance gain using
gsl_sf_pow_int() in place of pow(), even for very large powers.  GSL
developers might keep that in mind when performing x^n type
computations.  By the way, I noticed GSL also has a gsl_pow_int()
with essentially the same implementation in the sys/pow_int.c file;
which is preferred?

regards,
jt

[*] Kachitvichyanukul, V. and Schmeiser, B. W.  Binomial Random
   Variate Generation.  Communications of the ACM, 31, 2 (February,
   1988) 216. There is also an algorithm called PTPE for Poisson
   variates; Dr. Schmeiser gave me a fortran implementation and I
   have GSL'ified that as well, and will be posting it soon.

---------------------------------------------
James Theiler                     jt@lanl.gov
MS-B244, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

Attachment: new-binomial.diff
Description: Text document

Attachment: new-binomial.tar.gz
Description: application/gunzip