cosh magic number?
brde at optusnet.com.au
Fri May 31 20:51:34 UTC 2013
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
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.
More information about the freebsd-numerics