Pull request #33 on the `StatsFuns`

package provides a link to Catherine Loader’s `note`

on evaluating the binomial and Poisson probability mass function. One central part to this is evaluating the error in Stirling’s approximation to the log-factorial function. The version used is odd powers of `1/n`

with coefficients that are inverses of integers. It starts as

Are there clever ways of evaluating such a polynomial so that, e.g., switching to BigInt would be transparent? The code used by Catherine Loader and translated more-or-less literally in the PR memoizes the values for 1 to 15 (where the error is greatest) then switches to Horner evaluation of polynomials with explicitly evaluated Float64 coefficients. I think it would be better to retain the integer representation of the inverse of the coefficients.