Date:Tue, 11 Jun 2013 11:22:09 +1000 (EST)From:Bruce Evans <brde@optusnet.com.au>To:Bruce Evans <brde@optusnet.com.au>Cc:David Schultz <das@freebsd.org>, freebsd-numerics@freebsd.org, Steve Kargl <sgk@troutmask.apl.washington.edu>Subject:Re: Implementation for coshl.Message-ID:<20130611104252.G50308@besplex.bde.org>In-Reply-To:<20130611040059.A17928@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>

Next in thread | Previous in thread | Raw E-Mail | Index | Archive | Help

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

Want to link to this message? Use this URL: <http://docs.FreeBSD.org/cgi/mid.cgi?20130611104252.G50308>