Implementation for coshl.

[Bruce Evans]

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