Date: Sun, 10 May 2015 13:59:46 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: small cleanup patch for e_pow.c Message-ID: <20150510205946.GA85935@troutmask.apl.washington.edu> In-Reply-To: <20150511054753.L1656@besplex.bde.org> 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>
next in thread | previous in thread | raw e-mail | index | archive | help
On Mon, May 11, 2015 at 06:36:51AM +1000, Bruce Evans wrote: > On Sun, 10 May 2015, Steve Kargl wrote: > > > On Sun, May 10, 2015 at 08:16:14PM +1000, Bruce Evans wrote: > >> ... > >> So 1 is the only numbers near 1 that doesn't give overflow. > > > > Thanks for the explanation! That help dislodge a mental block. > > > > 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. (1-2**(-p))**(2**N) = 2**(emin-p) Using log() and log(1+z) = z for |z| << 1, the above equation gives N = p + log((p-emin) * log(2)) / log(2) where the equality is of course approximation. p = 24 --> N = 30.690 p = 23 --> N = 29.681 -- Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20150510205946.GA85935>