From owner-freebsd-numerics@FreeBSD.ORG Fri May 31 20:51:34 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.FreeBSD.org [8.8.178.115]) by hub.freebsd.org (Postfix) with ESMTP id A5DA542A for ; Fri, 31 May 2013 20:51:34 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail10.syd.optusnet.com.au (mail10.syd.optusnet.com.au [211.29.132.191]) by mx1.freebsd.org (Postfix) with ESMTP id 2EB08C01 for ; Fri, 31 May 2013 20:51:33 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail10.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id r4VKpOxf018305 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Sat, 1 Jun 2013 06:51:26 +1000 Date: Sat, 1 Jun 2013 06:51:24 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl Subject: Re: cosh magic number? In-Reply-To: <20130531191410.GA74343@troutmask.apl.washington.edu> Message-ID: <20130601052415.H15844@besplex.bde.org> References: <20130531191410.GA74343@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=e/de0tV/ c=1 sm=1 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=WUpannN5onsA:10 a=qPt0-ISivhtabmaQ0fEA:9 a=CjuIK1q_8ugA:10 a=gwKr3FwWfh0Jz3qo:21 a=etSL3OnsOLLAcWig:21 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 Cc: freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 31 May 2013 20:51:34 -0000 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 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