hello,

  This message is mostly for Gerard Jungman, since it relates to his
specfunc/gamma.c, but others may have worthwhile opinions as well.

  While looking at another part of GSL (randist/binomial.c -- I'll
have more to say about that in a later email), I realized that the
function gsl_sf_lnfact() [this computes log(n!)] could be speeded up
for small arguments by pre-computing the logarithms.  For small
arguments, the existing code already pre-computes real (double) and
integer (long) values, and maintains a static array of those values.
So including pre-computed logarithms is not a major design change.
For the current gsl_sf_lnfact() call, it takes the logarithm of the
precomputed value:
     result->val = log(fact_table[n].f);

  I altered the fact_table struct array to include an extra field,
double lnf; and then altered the array initialization to include the
logarithm of those first values.  So now the above line is replaced by
     result->val = fact_table[n].lnf;

  This is about 3x faster on my machine than taking the log.  (I was
hoping for about 12x improvement, based on initial experiments with a
fact_table array built directly into my test code, but there appears
to be a fair bit of overhead involved with going through the GSL -- I'm
guessing it has to do with the computation of the error estimates, which
are computed even if they are not used.  Error estimates are very worthwhile
things to have around, so we should think hard before we go in and try to
optimize them away.)

  While I was at it, I also modified the doub_fact_table in a similar way,
and the gsl_sf_lndoublefact() routine.  [This returns log(n!!)]

  I'll attach the diff file (I tried attaching the full gamma.c, but
the mail daemon bounced it because it was too big); I'll happily
send the full gamma.c to anybody who asks for it.  Meanwhile, I
have a few more general questions:

  1. I used a perl script to rebuild the precomputed fact_table[] and
     doub_fact_table[].  Should this perl script be part of GSL?  (Mark,
     you'll be relieved to know it's not part of the build script!)

  2. For the logarithms, we could easily extend the array well beyond the
     current bound of 170.  The cost is a large static array as part of GSL,
     and the question is: how large?  By the way, we could save some space
     by making the integer, double, and log tables separate.  One of length
     12, one of length 171, and one of substantially larger length?

  3. Speaking of saving space, do we really need the 'n' component of
     the fact_table struct?  It isn't used in gamma.c as far as I can
     tell.  To be sure, that makes it easier on the eye to cross-check
     with other tables that the values are correct, but we could just
     as well do something like
     static double fact_table_d[] = {
            /*  0 */  1.0,
            /*  1 */  1.0,
            /*  2 */  2.0,
            /*  3 */  6.0,
            ...etc
     }

  4. Also, a trick for computing the table size at compile time, so that
     you don't have to worry about simultaneously updating the
     '#define FACT_TABLE_MAX' and the table itself, is to use -- after
     the actual initialization such as that above for fact_table_d, a
     define statement along the lines of
        #define FACT_TABLE_D_SIZE (sizeof(fact_table_d)/sizeof(double))
        #define FACT_TABLE_D_MAX (FACT_TABLE_SIZE-1)
     Sometimes, I've also seen sizeof(fact_table_d[0]) in place of
     sizeof(double); one less place to hardcode the table's type.

regards,
jt

---------------------------------------------
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 -----