cosh magic number?

[Bruce Evans]

On Fri, 31 May 2013, Steve Kargl wrote:

> In msun/src/e_cosh.c, one finds the comment
>
> *
> *                                 exp(x) +  1/exp(x)
> * ln2/2 <= x <= 22 :  cosh(x) := -------------------
> *                                        2
>
> Where does the magic number 22 come from?

It is just a threshold at which a sloppier approximation becomes
adequate.  But you know that...

> Using exp(-|2x|) = 2**(1-p) with p = 53 for double, I
> arrive at 18.022, which is a little too small.

I get 18.368 using exp(-|2x|) = 2**p for the natural threshold.
(Consider x+y instead of E+1/E.  When x is 1+eps (with eps giving
1 in the last place, adding y = eps/2 causes rounding up to even).
This y is 2**p times smaller than x.  If x has extra precision,
then y still needs to start more than <full number of bits in x>
bits further out for adding y to have no effect, even if the
final result has no extra precision.)

I first thought that the extras are guard bits.  Perhaps they are,
but guard bits are not representable unless there is extra precision.
22 gives log2(exp(44)) ~= 64.479 bits.  64 fits well with x86 extra
precision.

This is easier to test in float precision.  Try all integer thresholds
near the chosen one, on all x.  Expect a difference for extra precision.

Bruce