Re: faster log factorials

[James Theiler <jt at lanl dot gov>]

On Mon, 12 May 2003, James Theiler wrote:

] On 12 May 2003, Gerard Jungman wrote:
]
] ] On Mon, 2003-05-12 at 08:02, James Theiler wrote:
] ]
] ]
] ] > A still-open question:  If we provide
] ] > pre-computed values, how many should we provide?  For the straight (no
] ] > logarithm) values, there is a natural cutoff at 170 since 170! (or is
] ] > it 171!?) is the largest value that is a valid IEEE double precison
] ] > number.  But for the logs, we can assume the higher the cutoff the
] ] > more often we'll be able to provide a fast precomputed value.  We
] ] > could easily provide thousands, and I think most computers nowadays
] ] > would not begrudge the memory.  But there may be other issues that I
] ] > am not considering.
] ]
] ] The natural cutoff is where Stirling's formula becomes good enough.
]
] good suggestion.
]
] ] Computing that way will also be better than a memory access,
] ] I imagine.
] ]
] ] If we start using Stirling at 170 then it should be good
] ] enough to take up to sixth order in the series correction.
] ] So we can probably use the same arrays, etc.
] ]
] ] I don't remember how the logic works in the factorial functions.
] ] It might need a little tweaking to get optimal performance,
] ] so it takes the right branch without too much fuss.
]
] For large arguments, seems to be gsl_sf_lnfact_e -> gsl_sf_lngamma_e
] -> lngamma_lanczos which implements a Stirling-like formula.
] However, Lanczos requires two log evals (cf Stirling requires one),
] and Lanczos and Stirling both require several terms of a polynomial
] (Lanczos has eight terms). Going straight to Stirling from the ln_fact
] code could save a few branches, and one logarithm, and possibly a few
] polynomial terms.  My intuition is that that won't beat a table
] lookup, but it could be a substantial improvement. (And we have to
] switch over from table lookup to large n formula at *some* point.) If
] I get a chance, I'll give it a try.
]
] jt

Hi Gerard and others,

I looked into timings for different ways of computing log(n!).

There are basically three approaches.  One is direct table lookup,
which obviously works only for n less than the size of the table. Two
is to call gsl_sf_lngamma(n+1.0), which currently employs a Lanczos
approximation for all arguments except n=0 and n=1.  A third approach,
suggested above, is to use the venerable Stirling's approximation [eg,
see Abramowitz and Stegen, 6.1.41].

Currently, gsl_sf_lnfact(n) uses table lookup for n<=170, and
the Lanczos formula for large n.  A relevant point for comparison,
therefore, is n ~ 170, and at that value I find for my machine [*]

Lookup:   0.05 usec/eval
Lanczos:  1.03 usec/eval
Stirling: 0.34 usec/eval

(usec = 1e-6 sec) Actually, the value for Stirling varies from
0.3-0.4, depending on the number of terms in the polynomial correction
are used.  I found that two terms were adequate for double precision
on my machine [*] for large numbers: n > 170.

That is, Lookup is 7x faster than Stirling, which in turn is 3x faster
than Lanczos.  To my mind, this is an argument in favor of running up
the upper bound on the table for log(n!); further, whenever n exceeds
whatever upper bound we choose, we should use Stirling instead of
Lanczos for the large n result.

In fact, we might even consider Stirling instead of Lanczos for large
arguments of the gsl_sf_lngamma() function as well.  Admittedly, this
does make the code that much more complicated...

regards,
jt

[*] my computer is a two-headed Intel Pentium linux box running at
933Mhz.  Only one head is used in these computations.  I would expect
the relative timings to be (more or less) invariant to computer type.

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