Date: Sat, 1 Jun 2013 06:51:24 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-numerics@FreeBSD.org Subject: Re: cosh magic number? Message-ID: <20130601052415.H15844@besplex.bde.org> In-Reply-To: <20130531191410.GA74343@troutmask.apl.washington.edu> References: <20130531191410.GA74343@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130601052415.H15844>