I don't know how many folks consider practical numerical computation to be "fun", but I do (-: I've been looking for a numeric Gamma method that gives 34 digits of precision suitable for quad-precision floating-point [2] although I would settle for 31 digits suitable for "double-double". On Mon, Jan 9, 2012 at 02:24, Bill Gosper <billgosper@gmail.com> wrote:
I missed where Warren's paper probably says that, for "unrestricted" (arbitrary precision) numerics, the divergence of the log Gamma and polygamma series is rather irrelevant. A plausible algorithm is simply to add some positive integer n to z so that the Stirling series z+n swoops within the desired accuracy. Then divide by (z+1)(z+2)...(z+n). Windschitl's and Warren's improvements make this even more attractive. [...]
In the 2006 paper from Warren Smith [1], (which is very good by the way, now that I've read it more carefully) I am only seeing Windschitl in one place, credited with the approximation: z! = z Gamma[z] = Sqrt[2*Pi*z] * [z/E]^z * [z*Sinh[1/z]]^[z/2] * [1+O(|z|^-5)] which (if the O(|z|^-5) part is taken at face value) means that for 34 digits, I'd need to ensure z>6300000 (6.3 million) and that clearly won't work. With the Stirling series for log(x!), the biggest term is x log(x), and the smallest is of magnitude BernoulliB[2*n]/[2*n*[2*n-1]*x^[2*n-1]] Back in 1998 I studied this for the case of 15 digits of precision, and I worked out that I could I sum 8 terms, with the first omitted term being the "+43867/244188*x^17" term, which gives 15 digits provided my argument x is at least 7. No problem, quite reasonable. But if I want near 34 digits of precision, it seems I need to sum the terms at least as far as the x^-23 term, with the residual being of order 657931/(300 x^25). With that number of terms, x needs to be at least 24. The worst case is when computing something like Gamma[0.6]. In that case I need to: - Multiply out the product P=(x+1)(x+2)... which has 23 terms (22 multiply plus 23 add) - Sum the terms of the series (12 terms: 12 add, 24 multiply and 12 divide) - compute Log[n] and Log[P] as always, plus the other lead terms (a few more primitive operations) That's an awful lot of operations, which is why I still want to use Lanczos' method instead, or Gosper's latest work if that turns out to give me the digits I need. The Lanczos method has the same requirement for two calls to Log[], but it appears I can get 34 digits of precision with about 17 terms, using 2 adds and one divide per term. For example, at [3] and [4] is a 15-digit precision version using 6 terms, which is 12 adds and 6 divides and Viktor Toth (see [5]) uses 80 terms to get about 165 digits of precision. I haven't multiplied out the matrices to find coefficients for 34-digit precision (not yet). Maybe someone can set me straight on this... - Robert More references and more flexible source code at: http://mrob.com/pub/ries/lanczos-gamma.html [1] http://rangevoting.org/WarrenSmithPages/homepage/gammprox.pdf [2] http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format [3] William Press et al., Numerical Recipes in C, section 6.1, function "gammln" (pp. 213-216 in the hardcover 2nd edition) [4] http://paulbourke.net/miscellaneous/functions/ (which gives the same code as [3]) [5] http://www.vttoth.com/CMS/projects/41-the-lanczos-approximation -- Robert Munafo -- mrob.com Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com