Skip site navigation (1)Skip section navigation (2)
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>