Le 05/02/2012 14:43, Dima Pasechnik a écrit :
On Sunday, 5 February 2012 21:03:46 UTC+8, Snark wrote:
(b) 1, 2 and 4 are all the same bug in various disguise : the libc
lgammal function isn't good enough ;
more precisely, the (e)glibc gamma function is not good. (lgammal isn't
too good either, but it's a different story).
Let me post your link to a good implementation of gamma(), in netbsd's libc:
http://cvsweb.netbsd.org/bsdweb.cgi/~checkout~/src/lib/libm/noieee_src/n_gamma.c?rev=1.7&content-type=text/plain&only_with_tag=MAIN
(and it states the error bounds, by the way!)
(I meant gammal and wrote lgammal indeed)
Yes, the code is wonderful, but isn't for double-precision I think. And
there's more ; since I was looking for a good and portable version, I
turned to NetBSD ; and since I didn't find exactly what I wanted, I
joined #netbsd to discuss, and I learned :
- 4.2BSD introduced the wrongly named gamma function as being the
logarithm of the gamma function ;
- in 4.3BSD, it got renamed to lgamma, but gamma was still an alias to it ;
- it's only in 4.4BSD that things were corrected, but NetBSD still uses
the 4.3BSD convention, even though they're supposed to be 4.4BSD, see
[1] for example (I'm leaving a trail of bug reports everywhere I set a
foot these days...)!
- even for the GNU libc, things are muddy, as you can see in [2], where
again they mention that gamma is a lgamma!
- in C99, the gamma function is supposed to be named tgamma if I
understand well, perhaps precisely because by the above points, the
situation is... interesting!
SSSSSSOOOOOOOOO... it seems trying to use "gamma" from a system
libc/libm is just asking for pain. Dima proposed on irc to just use gsl,
and I would agree:
(1) it's already there and
(2) it will shield sage nicely from what looks like a historical conundrum.
this is explicit in 1 and 2, and
implicit in 4 where for floats the computation is by a quotient of
gamma
calls [aside: that is probably a wrong way to do things].
I'm sure there is research done on how to do this numerically in a
proper way...
Well, even without going to the research level, using
gamma(x+1)=x*gamma(x) makes it easy to do computations for
gamma(a)/gamma(b) with numbers big like |a-b| instead of like max{a,b}.
But that's unrelated with that thread and I opened a ticket about it.
(c) 3 is a different problem, and although "#tol rel" makes it pass the
tests, the fact that it needs it to pass is bad.
So my proposition is the following:
(1) I'll try to fix the implementation of lgammal upstream ;
gammal, and gamma, the first thing, perhaps. Not sure if lgammal can be
fixed easily, actually.
See above : there's courage and there's inconscience. Let's flee!
(2) I'll try to investigate the maxima string-to-float conversion ;
(3) the sage developpers should still decide and make public some
statement as to what they consider a correct numerical approximation,
because that thread only made obvious that there's a need for a
reference document.
This is a very good idea. The link to netbsd's implementation show
that's it's doable...
Yes, that code is beautiful, its documentation is beautiful.
Snark
PS:
[1]
http://cvsweb.netbsd.org/bsdweb.cgi/~checkout~/src/lib/libm/src/w_gamma.c?rev=1.11&content-type=text/plain&only_with_tag=MAIN
[2]
http://www.gnu.org/software/libc/manual/html_node/Special-Functions.html#Special-Functions
--
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
