From owner-freebsd-numerics@FreeBSD.ORG Tue Jun 11 01:22:19 2013 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id 5299B413; Tue, 11 Jun 2013 01:22:19 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail108.syd.optusnet.com.au (mail108.syd.optusnet.com.au [211.29.132.59]) by mx1.freebsd.org (Postfix) with ESMTP id DCA5C1ED8; Tue, 11 Jun 2013 01:22:18 +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 mail108.syd.optusnet.com.au (Postfix) with ESMTPS id B175D1A06F7; Tue, 11 Jun 2013 11:22:10 +1000 (EST) Date: Tue, 11 Jun 2013 11:22:09 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Bruce Evans Subject: Re: Implementation for coshl. In-Reply-To: <20130611040059.A17928@besplex.bde.org> Message-ID: <20130611104252.G50308@besplex.bde.org> References: <20130610003645.GA16444@troutmask.apl.washington.edu> <20130610110740.V24058@besplex.bde.org> <20130610055834.GA68643@zim.MIT.EDU> <20130610162927.GA20856@troutmask.apl.washington.edu> <20130610170603.GA72597@zim.MIT.EDU> <20130611040059.A17928@besplex.bde.org> 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=K8x6hFqI c=1 sm=1 a=LM0AswAWfpYA:10 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=MtpAjnDH-vAA:10 a=Udj031w0ckm0zOnfpgkA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 Cc: David Schultz , freebsd-numerics@freebsd.org, Steve Kargl 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: Tue, 11 Jun 2013 01:22:19 -0000 On Tue, 11 Jun 2013, Bruce Evans wrote: > On Mon, 10 Jun 2013, David Schultz wrote: > >> On Mon, Jun 10, 2013, Steve Kargl wrote: >... >>> I originally had >>> >>> if (ix < BIAS - 33) { /* |x| < 0x1p-33 */ >>> /* cosh(+-0) = 1. */ >>> if (x == 0) >>> return (one); >>> /* cosh(tiny) = one */ >>> return (one + tiny * tiny); >>> } >>> >>> because C99 requires cosh(+-0) to be exact. For |x| != the >>> return value will raise inexact. > ... > No, the above is quite broken. tiny * tiny gives underflow, and > when the underflow is to 0, adding 1 to it doesn't give inexact. > Also, the 'tiny' in the comment is confusing, since it refers to > a different variable than the 'tiny' in the code. > > Just copy the known-good fdlibm code. It uses one+t, where t is > expm1(fabs(x)). This has many subtleties: > - it avoids pessimizing the usual case > - it moves all the complications for setting inexact to expm1(). These > include: > - the expression 1+x doesn't work for tiny x since it rounds incorrectly > in some modes. It also requires x >= 0 for correct rounding, so > we may have to caclulate |x| earlier than we want. > - the expression 1+x*x doesn't work for tiny x since it may underflow > (just like 1+tiny*tiny) > - the expression x doesn't work for tiny x in expm1() since it doesn't > set inexact. > - the expression x+x*x doesn't work for tiny x in expm1() since it may > underflow. Bah, this is more magic and broken than first appeared. The fdlibm code is now known-bad :-). I stumbled across this when implementing the use of k_expf() in ccoshf(). o One subtlety in the fdlibm cosh() is that it uses an artificially small threshold of 2**-55 to "ensure" that 1+this gives 1 with inexact (actually 1+eps in some rounding modes). A threshold of about 2**26.5 is sufficient for cosh() to decay to 1, but fdlibm wants to abuse the tiny value given by each arg below the threshold for setting inexact as a side effect of evaluating the expression 'one + expm1f(x)'. This obviously requires a threshold of <= 2**-53 or 2**-54. The value of 2**-55 actually used makes no sense. o I broke this in coshf() in 2005 when fixing its threshold. The threshold of 2**-55 was blindly copied. It should have been reduced to 2**-26 for a direct translation. I changed it to the natural threshold of 2**-12. This also gives a micro-optimization. This made extra-precision bugs more apparent, and I fixed them by returning 'one' instead of 'one + expm1f(x)' for tiny x. Returning 'one' gives another micro-optimization. But this weakens the support for non-default rounding modes. o fdlibm cosh()'s large threshold is not large enough to actually work. 'one + 2**-55' may be evaluated and returned in extra precision, and it only takes 2 or 3 bits extra for it to evaluate to itself instead of 1. coshf()'s blindly copied threshold accidentally protected against using '1.0F + 2**-26F' which always evaluates to itself on i386, and is then normally returned unchanged. Assertions like assert(cosh(tiny) == 1) should then fail, but cosh() may return extra precision so is also a bug in the assertions -- essentially the same one that I found recently for exp() non-overlow. This is still a bug in cosh() which would be found by a more careful assertion -- cosh(x) may be evaluated in extra precision, but then its value must be ~1+x*x/2, not ~1+|x|. FreeBSD's default rounding precision protects against these bugs in both float and double precision provided the theshold is <= 2**-53 for both. Bruce