I think we may be overlooking a very reasonable option. Python already has a gamma function in the math module! It is a separate implementation that does not depend on libc, and it gives reasonable results (though perhaps not as good as eglibc tgammal() on x86):
sage: max(( abs( RR(math.gamma(float(x))) - x.gamma())/x.gamma().ulp() for x in (RR.random_element(-171, 171) for _ in xrange(100000)) )) 6.00000000000000 sage: mean([ abs( RR(math.gamma(float(x))) - x.gamma())/x.gamma().ulp() for x in (RR.random_element(-171, 171) for _ in xrange(100000)) ]) 0.867960000000000 sage: median([ abs( RR(math.gamma(float(x))) - x.gamma())/x.gamma().ulp() for x in (RR.random_element(-171, 171) for _ in xrange(100000)) ]) 1.00000000000000 sage: sage: [ZZ(math.gamma(float(n))) == n.gamma() for n in srange(1, 23)] == [True] * 22 True The source code does say: In extensive but non-exhaustive random tests, this function proved accurate to within <= 10 ulps across the entire float domain. Note that accuracy may depend on the quality of the system math functions, the pow function in particular. So if the accuracy of pow() in eglicb relies on long doubles, then there may be a problem, but maybe it will work well there. On Sun, Feb 5, 2012 at 2:16 PM, Jonathan Bober <jwbo...@gmail.com> wrote: > See http://trac.sagemath.org/sage_trac/ticket/12449 > > I'm in the middle of [hopefully] fixing this by calling the gsl gamma > function. While I'm at it, I'll also make the evaluation on basic types > much faster, as it shouldn't go though Ginac. (Actually, I've already > mostly written a fix. I opened the ticket and wrote this email while > running 'make ptestlong'.) > > Have you actually checked the gsl implementation on ARM? For me it at > least satisfies > > sage: [ZZ(sage.gsl.special_functions.gamma(n)) == n.gamma() for n in > srange(1, 23)] == [True] * 22 > True > > (Unfortunately, I just wrote that sage.gsl.special_functions.gamma(). > There is no high level interface for you to immediately check if it gives > the right answer.) > > ----- > > Never mind all that: the gsl implementation is not very good at all, > whereas the libc implementation on my machine seems quite good. Old (libc): > > sage: gamma(1.23).str(base=2) > '0.11101001001001110011101011110010110010101100101100001' > sage: RR(gamma(1.23r)).str(base=2) > '0.11101001001001110011101011110010110010101100101100001' > > gsl: > > sage: RR(gamma(1.23r)).str(base=2) > '0.11101001001001110011101011100000000000111110110111001' > > look at all the wrong digits! In my testing, gsl_gamma() has a typical > error of around 1e-8 on random input between 1 and 2, the tgammal() rounded > to double precision has a typical error of 0 (compared to the correctly > rounded value from mpfr). > > What do you get on ARM if you do something like > > sage: for n in range(100): > ....: x = RR.random_element(1, 2) > ....: print abs(RR(gamma(float(x))) - x.gamma()) > ....: > > ? > -- To post to this group, send an email to sage-devel@googlegroups.com To unsubscribe from this group, send an email to sage-devel+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URL: http://www.sagemath.org
