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

["Joseph S. Myers" <joseph at codesourcery dot com>]

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.

> > +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).

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

-- 
Joseph S. Myers
joseph@codesourcery.com