On Thursday, 2 February 2012 12:52:53 UTC+8, Jonathan Bober wrote:
>
> 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. */
> [...]
>
>
yes, that's correct, I think (at least I checked this to be the case on
ARM).
> 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.
>
here, for x=6, the relative error is >10^(-15).
for x=150, it gets bigger, > 10^(-13).
> (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 <> 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
