Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)

["Carlos O'Donell" <carlos at redhat dot com>]

On 05/03/2013 06:30 AM, Joseph S. Myers wrote:
> On Thu, 2 May 2013, Carlos O'Donell wrote:
> 
>> How did you choose the number of digits in the answer?
> 
> I generally use 40 decimal digits, which is safely more than enough for 
> any supported floating-point format.
> 
>> How did you compute the answer?
> 
> Using MPFR.
> 
>> What am I misunderstanding about the way libm-test.inc works?
> 
> In general, the expected results given are decimal versions of the 
> mathematically exact result, rounded to some number of decimal digits.  
> The number should be sufficient that when GCC then rounds to the 
> particular floating-point format being tested, the double rounding is 
> unlikely to have results different from single rounding direct to the 
> appropriate (binary) format.

That is as I understood it, but it would seem to me that 40 digits 
isn't enough. Did I do something wrong?

>>> +static const long double gamma_coeff[] =
>>> +  {
>>> +    0x1.5555555555555555555555555555p-4L,
>>> +    -0xb.60b60b60b60b60b60b60b60b60b8p-12L,
>>> +    0x3.4034034034034034034034034034p-12L,
>>> +    -0x2.7027027027027027027027027028p-12L,
>>> +    0x3.72a3c5631fe46ae1d4e700dca8f2p-12L,
>>> +    -0x7.daac36664f1f207daac36664f1f4p-12L,
>>> +    0x1.a41a41a41a41a41a41a41a41a41ap-8L,
>>> +    -0x7.90a1b2c3d4e5f708192a3b4c5d7p-8L,
>>> +    0x2.dfd2c703c0cfff430edfd2c703cp-4L,
>>> +    -0x1.6476701181f39edbdb9ce625987dp+0L,
>>> +    0xd.672219167002d3a7a9c886459cp+0L,
>>> +    -0x9.cd9292e6660d55b3f712eb9e07c8p+4L,
>>> +    0x8.911a740da740da740da740da741p+8L,
>>> +    -0x8.d0cc570e255bf59ff6eec24b49p+12L,
>>> +  };
>>
>>
>> If all three of these tables are of type long double, why can't
>> we define the table once and use it everywhere?
> 
> Each table has as contents numbers with the appropriate number of hex 
> digits for the applicable long double format (rounded once by MPFR doing a 
> conversion from an exact rational (GMP mpq)).  And the length of each 
> table, and the associated threshold above which Stirling's approximation 
> is used, were chosen separately for each format.  Although the first issue 
> could also be handled by having precise-enough values to round correctly 
> for every format (like in libm-test.inc), a single copy would be more 
> awkward for varying the length of the table in each case.  Only one copy 
> of the table is used in any particular build of glibc.  And this table is, 
> in any case, small (compared for example to the 35 degree-19 polynomials I 
> referred to used by the IA64 lgammal to expand about the roots of lgamms).

Thanks, that makes sense.

> Recall that the series in Stirling's approximation is divergent 
> everywhere, but the error when it is truncated has the same sign, and 
> lesser magnitude, than the first neglected term.  To be able to use this 
> approximation for a given format and input value, the series terms must 
> drop small enough for that input value before they start to rise.  This 
> places an absolute lower bound on where the approximation can be used for 
> a given floating-point format - and then it makes sense in fact to start 
> using the series only at a higher value, trading off using fewer terms of 
> the series with more uses of the recurrence being needed to get into the 
> Stirling range.  (E.g., for ldbl-128 it's possible to use the series for 
> arguments at least 13, but then you need to use terms up to the 53rd 
> power, whereas by starting using it at 24 only terms up to the 27th power 
> are needed.)

That makes sense, and limits the table size correct?

Cheers,
Carlos.