Re: [math-fun] Gamma function (formulas collection)
(Sorry about (no subject).) Hmm, "only" 26 digits with seven b*ggerfactors on z>9: http://gosper.org/munafacmu26.png Note that the z^7 means we only need ~20 digit values for the seven b*ggerfactors, but it looks like we'll need > eight, ior larger "9" for quad (34digit) precision. How many divides is a "+a5*z^2" worth? --rwg On Tue, Jan 10, 2012 at 4:12 AM, Bill Gosper <billgosper@gmail.com> wrote:
Robert>
I only need it to work for real arguments. But if that were an issue I could stick with Lanczos.
Any given set of Lanczos parameters work equally well on the entire half-plane Real[z]>0, and with the same guaranteed error bound. This article describes it really well:
http://www.haoli.org/nr/bookcpdf/c6-1.pdf
which is from: William H. Press et al., Numerical recipes in C: the art of scientific computing (2nd edition). Cambridge University Press, (1992) ISBN 0521431085.
On Mon, Jan 9, 2012 at 12:06, Warren Smith <warren.wds@gmail.com <http://gosper.org/webmail/src/compose.php?send_to=warren.wds%40gmail.com>> wrote:
Do you want it to work for complex or only for real argument?>> > 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".> -- 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
Just two b*ggerfactors, da*n near double precision on z>9: http://gosper.org/munafacmu16.png Probably could be double prec for z>10. (Manually Remezed coefficients aren't finished, but mu won't change much.)
I'm guessing quad will take eight b*ggerfactors. --rwg
I have gotten the 10-term "munafacmu34" version performing to near the limits of what's possible in a 32-digit floating-point (Dekker style double-double) numerics library. To answer your question *"How many divides is a '+a5*z^2' worth?"*, I'd say it's probably about a one-to-one tradeoff. In other words, if you add one term to a polynomial to save one divide in the "Gamma[x+1]=x Gamma[x]" iteration, that would be about a break-even. So if you can reduce the minimum argument by 2 via adding one term, that's a win. The bigger issue I'm dealing with is the transcendental operations. My Sinh() is doing an 8-step Taylor series (for 1/z <= 1/14), Exp() is a 16 to 20 step Taylor series, and Log() is a little slower than Exp(). Since we need about 5 of these transcendentals to use the approximation formula, that's where most the computation time is being spent. - Robert [1] See bibliography at mrob.com/pub/math/f161.html On Tue, Jan 10, 2012 at 19:14, Bill Gosper <billgosper@gmail.com> wrote:
(Sorry about (no subject).) Hmm, "only" 26 digits with seven b*ggerfactors on z>9:
http://gosper.org/munafacmu26.png
Note that the z^7 means we only need ~20 digit values for the seven b*ggerfactors, but it looks like we'll need > eight, ior larger "9" for quad (34digit) precision. How many divides is a "+a5*z^2" worth? --rwg
-- 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
participants (2)
-
Bill Gosper -
Robert Munafo