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

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

On Fri, 3 May 2013, Carlos O'Donell wrote:

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

For ldbl-96, 21 digits is enough if the decimal value is the 
correctly-rounded decimal version of a representable floating-point number 
(that's the value of DECIMAL_DIG).  You appear to be giving the 
full-precision decimal version of a particular floating-point number, but 
that's not necessary for the compiler to round correctly to a uniquely 
determined floating-point number.

Of course, since the number given here isn't exactly representable - it's 
a rounded version of the actual value of the gamma function, not of a 
64-bit approximation to the actual value of the gamma function - there is 
no a priori bound on how many digits is enough, since a value could be 
arbitrarily close to half way between two representable values.  (strtod 
is only able to stop after a certain point based on knowledge of whether 
any of the digits that follow are 0.)  For measuring accuracy of results 
of most libm functions, a maybe 1 in 100000 chance for ldbl-128, much less 
for other formats, that the decimal value rounds wrongly and so the 
reported ulps differ from the precise ulps by plus or minus 1 - in a case 
where the actual number of ulps, relative to the exact mathematical result 
rather than a rounded version thereof, is very close to half way between 
the two possible integer numbers of ulps - is not considered significant.

(Cf. <http://sourceware.org/ml/libc-alpha/2012-05/msg02040.html>, where 
some constants in <math.h> needed up to 36 digits after the point to get 
correct roundings for ldbl-128.)

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

Yes.

-- 
Joseph S. Myers
joseph@codesourcery.com