From owner-freebsd-numerics@FreeBSD.ORG Sun May 10 20:59:48 2015 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id 7BD238D6 for ; Sun, 10 May 2015 20:59:48 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 5AF901A4B for ; Sun, 10 May 2015 20:59:48 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id t4AKxkpl085971 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Sun, 10 May 2015 13:59:46 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id t4AKxkQS085970; Sun, 10 May 2015 13:59:46 -0700 (PDT) (envelope-from sgk) Date: Sun, 10 May 2015 13:59:46 -0700 From: Steve Kargl To: Bruce Evans Cc: freebsd-numerics@freebsd.org Subject: Re: small cleanup patch for e_pow.c Message-ID: <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> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20150511054753.L1656@besplex.bde.org> User-Agent: Mutt/1.5.23 (2014-03-12) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.20 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: Sun, 10 May 2015 20:59:48 -0000 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