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.
(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.)
On Wed, Feb 1, 2012 at 8:04 PM, Dima Pasechnik <dimp...@gmail.com> wrote:
>
>
> On Thursday, 2 February 2012 06:24:18 UTC+8, Robert Bradshaw wrote:
>
>> On Wed, Feb 1, 2012 at 4:19 AM, Julien Puydt <> wrote:
>> > ============= Forewords =============
>> >
>> > I investigated the numerical issues on my ARM build, and after much
>> poking
>> > around and searching, I found that I was chasing the dahu : the tests
>> were
>> > wrong, and the result were good.
>>
>> No, the tests are fine, the results are clearly wrong.
>>
>> > Let's consider the numerical failures one by one :
>>
>> The following comments apply to all 4 tests.
>>
>> > ============= 1/4 =============
>> > File "/home/jpuydt/sage-4.8/devel/**sage-main/sage/functions/**other.py",
>> line
>> > 511:
>> > sage: gamma1(float(6))
>> > Expected:
>> > 120.0
>> > Got:
>> > 119.99999999999997
>> >
>> > Let's see how bad this is :
>> > sage: res=gamma(float(6))
>> > sage: res
>> > 119.99999999999997
>> > sage: n(res,prec=57)
>> > 120.0000000000000
>> > sage: n(res,prec=58)
>> > 119.99999999999997
>>
>> sage: a = float(120) - 2^-46; a, a == 120
>> (119.99999999999999, False)
>> sage: a = float(120) - 2^-46; a, a == 120
>> (119.99999999999999, False)
>> sage: n(a, prec=57)
>> 120.0000000000000
>> sage: n(a, prec=57) == 120
>> False
>> sage: n(a, prec=57).str(truncate=False)
>> '119.9999999999999858'
>>
>> The string representation of elements of RR truncates the last couple
>> of digits to avoid confusion, but floats do not do the same. Changing
>> tests like this to have n(..., prec=...) and relying on string
>> formatting only confuses the matter (as well as making things ugly).
>>
>> It looks like ARM's libc returns incorrect (2 ulp off) values for
>> gamma(n) for small n (at least). This should be fixed, not hidden, or
>> at least marked as a known bug on ARM (only). IMHO this error is a
>> much bigger issue than the noise due to a numerical integration
>> arising from double-rounding in moving floats in and out of 80-bit
>> registers and other low-level details. That's what the tolerance
>> makers should be used for.
>>
>> > ============= CONCLUSION =============
>> > A double precision floating point number is supposed to have 53 digits,
>> > according to the norm
>> > (http://en.wikipedia.org/wiki/**IEEE_754-2008<http://en.wikipedia.org/wiki/IEEE_754-2008>),
>> and the
>> > results are correct from that point of view.
>>
>> No, their last (two) binary digits are wrong. If the test was
>>
>
> further digging shows that it's the implementation of gammal() in the
> platform's libc (eglibc is used)
> is to blame; they do expl(lgammal()), leading to loss of precision, as
> platform's long double is only 8 bytes, and it's simply not
> possible to stuff enough precision in log(gamma) if you only have 8 bytes!
>
> Dima
>
>
> sage: gamma(float(6)) == 120
>> True
>>
>> It would fail just as well.
>>
>> > So the tests should be modified not to depend on the specific
>> implementation
>> > : they're currently testing equality of floats!
>> >
>> > I would provide a patch for the tests so they use n(..., prec=53), but
>> I'm
>> > hitting a problem in one of the cases ; see the mail I sent yesterday
>> for
>> > more about that.
>>
>> See above.
>>
>> - Robert
>>
>> --
> 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
>
--
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
