This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: Carlos O'Donell <carlos at redhat dot com>
- Cc: <libc-alpha at sourceware dot org>
- Date: Fri, 3 May 2013 14:07:26 +0000
- Subject: Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)
- References: <Pine dot LNX dot 4 dot 64 dot 1305022233490 dot 12072 at digraph dot polyomino dot org dot uk> <51831AE0 dot 5020901 at redhat dot com> <Pine dot LNX dot 4 dot 64 dot 1305031016200 dot 26911 at digraph dot polyomino dot org dot uk> <5183A73F dot 70209 at redhat 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
- References:
- Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)
- Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)
- Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)
- Re: Improve tgamma accuracy (bugs 2546, 2560, 5159, 15426)