Date: Tue, 11 Jun 2013 11:22:09 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Bruce Evans <brde@optusnet.com.au> Cc: David Schultz <das@freebsd.org>, freebsd-numerics@freebsd.org, Steve Kargl <sgk@troutmask.apl.washington.edu> Subject: Re: Implementation for coshl. Message-ID: <20130611104252.G50308@besplex.bde.org> In-Reply-To: <20130611040059.A17928@besplex.bde.org> References: <20130610003645.GA16444@troutmask.apl.washington.edu> <20130610110740.V24058@besplex.bde.org> <20130610055834.GA68643@zim.MIT.EDU> <20130610162927.GA20856@troutmask.apl.washington.edu> <20130610170603.GA72597@zim.MIT.EDU> <20130611040059.A17928@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, 11 Jun 2013, Bruce Evans wrote:
> On Mon, 10 Jun 2013, David Schultz wrote:
>
>> On Mon, Jun 10, 2013, Steve Kargl wrote:
>...
>>> I originally had
>>>
>>> if (ix < BIAS - 33) { /* |x| < 0x1p-33 */
>>> /* cosh(+-0) = 1. */
>>> if (x == 0)
>>> return (one);
>>> /* cosh(tiny) = one */
>>> return (one + tiny * tiny);
>>> }
>>>
>>> because C99 requires cosh(+-0) to be exact. For |x| != the
>>> return value will raise inexact.
> ...
> No, the above is quite broken. tiny * tiny gives underflow, and
> when the underflow is to 0, adding 1 to it doesn't give inexact.
> Also, the 'tiny' in the comment is confusing, since it refers to
> a different variable than the 'tiny' in the code.
>
> Just copy the known-good fdlibm code. It uses one+t, where t is
> expm1(fabs(x)). This has many subtleties:
> - it avoids pessimizing the usual case
> - it moves all the complications for setting inexact to expm1(). These
> include:
> - the expression 1+x doesn't work for tiny x since it rounds incorrectly
> in some modes. It also requires x >= 0 for correct rounding, so
> we may have to caclulate |x| earlier than we want.
> - the expression 1+x*x doesn't work for tiny x since it may underflow
> (just like 1+tiny*tiny)
> - the expression x doesn't work for tiny x in expm1() since it doesn't
> set inexact.
> - the expression x+x*x doesn't work for tiny x in expm1() since it may
> underflow.
Bah, this is more magic and broken than first appeared. The fdlibm code
is now known-bad :-). I stumbled across this when implementing the use of
k_expf() in ccoshf().
o One subtlety in the fdlibm cosh() is that it uses an artificially small
threshold of 2**-55 to "ensure" that 1+this gives 1 with inexact
(actually 1+eps in some rounding modes). A threshold of about 2**26.5
is sufficient for cosh() to decay to 1, but fdlibm wants to abuse the
tiny value given by each arg below the threshold for setting inexact
as a side effect of evaluating the expression 'one + expm1f(x)'. This
obviously requires a threshold of <= 2**-53 or 2**-54. The value of
2**-55 actually used makes no sense.
o I broke this in coshf() in 2005 when fixing its threshold. The
threshold of 2**-55 was blindly copied. It should have been reduced
to 2**-26 for a direct translation. I changed it to the natural
threshold of 2**-12. This also gives a micro-optimization. This
made extra-precision bugs more apparent, and I fixed them by
returning 'one' instead of 'one + expm1f(x)' for tiny x. Returning
'one' gives another micro-optimization. But this weakens the support
for non-default rounding modes.
o fdlibm cosh()'s large threshold is not large enough to actually work.
'one + 2**-55' may be evaluated and returned in extra precision, and
it only takes 2 or 3 bits extra for it to evaluate to itself instead
of 1. coshf()'s blindly copied threshold accidentally protected
against using '1.0F + 2**-26F' which always evaluates to itself on
i386, and is then normally returned unchanged. Assertions like
assert(cosh(tiny) == 1) should then fail, but cosh() may return
extra precision so is also a bug in the assertions -- essentially
the same one that I found recently for exp() non-overlow. This
is still a bug in cosh() which would be found by a more careful
assertion -- cosh(x) may be evaluated in extra precision, but then
its value must be ~1+x*x/2, not ~1+|x|.
FreeBSD's default rounding precision protects against these bugs in both
float and double precision provided the theshold is <= 2**-53 for both.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130611104252.G50308>