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

Reply via email to