Le 02/02/2012 05:52, Jonathan Bober a écrit :
I've just been looking at this trying to figure out what was going on
and I was just going to say exactly the same thing.

I don't really know anything about the whole glibc vs eglibc thing, but
I bet the implementation is the same as
glibc-2.14.1/sysdeps/ieee754/dbl-64/e_gamma_r.c:

double
__ieee754_gamma_r (double x, int *signgamp)
{
   /* We don't have a real gamma implementation now.  We'll use lgamma
      and the exp function.  But due to the required boundary
      conditions we must check some values separately.  */
[...]

And sure enough:

sage: exp(RR(log(120))).str(truncate=False)
'119.99999999999997'

I'm not completely convinced that the results are wrong. They are
certainly less precise than the correct answers, and in this case the
correct answers can be represented exactly in double precision, but I
would not be surprised if on x86_64 there are still cases where the
returned answer is not as exact as possible. And as far as I can tell,
the underlying tgamma() works as advertised, in the sense that I cannot
find any description at all of how accurate it's answers are supposed to
be, so within 2ulp seems possibly reasonable to me.

Let us consider :
#include <math.h>
#include <stdio.h>

int
main ()
{
  long double x;
  x=6.0;
  printf ("%.20Lf\n", tgammal(x));
  x=10.0;
  printf ("%.20Lf\n", tgammal(x));
  return 0;
}

On an x86_64 box, I get :
$ ./test
119.99999999999999998612
362880.00000000000022737368

and on the ARM box, I get:
$ ./test
119.99999999999997157829
362880.00000000046566128731

Notice that even for the computation of tgammal(10), as I mentioned yesterday, the relative error is 1.2832....e-15, so...

(The maxima conversion is a different issue. There the result _is_
wrong. 1.7e17 is exact in double precision, and the conversion should
not change that.)

Well, that problem goes away with "# tol rel", hence the situation isn't that clear-cut.

Snark on #sagemath

--
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