Date: Mon, 11 May 2015 09:09:46 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-numerics@freebsd.org Subject: Re: small cleanup patch for e_pow.c Message-ID: <20150511084012.R2314@besplex.bde.org> In-Reply-To: <20150510205946.GA85935@troutmask.apl.washington.edu> References: <20150510002910.GA82261@troutmask.apl.washington.edu> <20150510113454.O841@besplex.bde.org> <20150510061458.GA82518@troutmask.apl.washington.edu> <20150510172810.E1812@besplex.bde.org> <20150510190114.GA85376@troutmask.apl.washington.edu> <20150511054753.L1656@besplex.bde.org> <20150510205946.GA85935@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Sun, 10 May 2015, Steve Kargl wrote: > On Mon, May 11, 2015 at 06:36:51AM +1000, Bruce Evans wrote: >> On Sun, 10 May 2015, Steve Kargl wrote: >>> >>> To find the magic numbers, it seems I need to consider >>> >>> (1-2**(-p))**(2**N) = 2**(emin-p) for underflow >>> (1+2**(-p))**(2**N) = (1-2**(-p))*2**emax for overflow >>> >>> With p = [24, 53, 64, 113], emin = [-125, -1021, -16381, -16381], >>> emax = [128, 1024, 16384, 16384], and the use of log(1+z) = z for >>> |z| << 1, I find >>> >>> underflow: N = [30.7, 62.5, 77.5, 126.5] >>> overflow: N = [30.5, 62.5, 77.5, 126.5] >> >> I plugged some numbers into pari to get sloppy estimates. E.g., >> (1 + 2^-23.0)^(2^N) / 2^128; then bump up N until the result is > 1. >> I get this at N = 29.5, so 30.7 seems too high. > > I used p = 24 in my (1+-2^(-p)). It seems that you're using p-1=23. It should be 1+2**-(p-1) and 1-2**-p, but the 1-epsilon side dominates; also 2**128 is sloppy with gradual underflow; in the scale factor for the 1-epsilon side, 128 must be replaced by something like 150. So 30.7 is correct if we don't bother to treat the sides differently. There is a related problem with extra precision. Expressions like huge*huge and tiny*tiny are always wrong in the presence of extra precision/exponent range, since they don't overflow or underflow but just give a wrong result if extra exponent range actually works. But extra precision/exponent is often broken: - C11 broke it and also broke efficiency by requiring extra precision/ exponent range to be destroyed on function return. clang doesn't implement this breakage, but gcc-4.8 does (with -std=c99 but not with -std=gnu99) - at least old versions of clang miscompile expressions like huge*huge unless huge is declared with the volatile hack. The C11 breakage makes returning the value INFINITY correct and then the only error in evaluating this at compile time is not raising the overflow exception. - at least old versons of both clang and gcc miscompile expressions like tiny*tiny, except... Functions like exp() are especially affected by this bug. cexp() is most noticeably broken by it, since it uses expressions like exp(x) * cos(y), written otherwise carefully by scaling exp(x) so as to avoid spurious overflow, so some non-overflowing cases return garbage related to huge*huge in extra precision. We mostly don't worry about this bug, except I fixed it for cexp*(). My test don't see it, since they assign function results to a variable without the extra precision, so as to destroy it like the C11 breakage requires. The bug should be fixed someday by replacing code like 'return huge*huge;' with 'return Inf_with_OE();'. Cases where the expression is actually x*x where x is known to be >= huge are harder to fix. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20150511084012.R2314>