From owner-freebsd-numerics@freebsd.org Thu May 18 23:47:30 2017 Return-Path: Delivered-To: freebsd-numerics@mailman.ysv.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by mailman.ysv.freebsd.org (Postfix) with ESMTP id 2DD05D73CB5 for ; Thu, 18 May 2017 23:47:30 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) 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 14C831763 for ; Thu, 18 May 2017 23:47:30 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.15.2/8.15.2) with ESMTPS id v4INlMhu077565 (version=TLSv1.2 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Thu, 18 May 2017 16:47:22 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.15.2/8.15.2/Submit) id v4INlMRp077564; Thu, 18 May 2017 16:47:22 -0700 (PDT) (envelope-from sgk) Date: Thu, 18 May 2017 16:47:22 -0700 From: Steve Kargl To: Bruce Evans Cc: freebsd-numerics@freebsd.org Subject: Re: Implementation of half-cycle trignometric functions Message-ID: <20170518234722.GB77471@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20170427231411.GA11346@troutmask.apl.washington.edu> <20170428175851.A1386@besplex.bde.org> <20170516224618.GA40855@troutmask.apl.washington.edu> <20170517094848.A52247@besplex.bde.org> <20170517180927.GA54431@troutmask.apl.washington.edu> <20170518072636.E951@besplex.bde.org> <20170518014226.GA63819@troutmask.apl.washington.edu> <20170518154820.A8280@besplex.bde.org> <20170518175434.GA74453@troutmask.apl.washington.edu> <20170519074512.M884@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20170519074512.M884@besplex.bde.org> User-Agent: Mutt/1.7.2 (2016-11-26) X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.23 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: Thu, 18 May 2017 23:47:30 -0000 On Fri, May 19, 2017 at 08:56:27AM +1000, Bruce Evans wrote: > On Thu, 18 May 2017, Steve Kargl wrote: > > > On Thu, May 18, 2017 at 04:34:57PM +1000, Bruce Evans wrote: > >> On Wed, 17 May 2017, Steve Kargl wrote: > >> > >>> On Thu, May 18, 2017 at 09:25:58AM +1000, Bruce Evans wrote: > >>>> On Wed, 17 May 2017, Steve Kargl wrote: > >>> ... > >>>>> As such, I've added a comment > >>>>> > >>>>> /* > >>>>> * pi_hi contains the leading 56 bits of a 169 bit approximation for pi. > >>>>> */ > >>>> > >>>> Why 56? 53 + 56 = 109. > >>> > >>> This is ld128. > >> > >> ld128 has precision 113. > > > > Yes. I know. > > > >>> static const long double > >>> pi_hi = 3.14159265358979322702026593105983920e+00L, > >>> pi_lo = 1.14423774522196636802434264184180742e-17L; > >>> > >>> pi_hi has 113/2 = 56 bits. > >>> pi_lo has 113 bit. > >>> > >>> 56 + 113 = 169 > >> > >> But hi is set intentionally sloppily by converting to double, so it has > >> 53 bits instead of the 57 that 56 would go with (57 = 113/2 rounded up). > > > > A 53-bit entity times a 56-bit entity needs 109-bits for the result. > > A 113-bit entity can hold a 109-bit entity. Or am I missing something? > > It wastes 4 bits. This reader has to guess if this is intentional for > magic elsewhere, or just sloppy. > > Also, the rounded value might naturally have some extra low zero bits. > So looking at the bits, it is not easy to see if there are the right > number. This is not worth a comment. The reader should trust that > the right (maximal) number of bits are used, and only check it > approximately. Unfortunately, I don't have the capibility to easily convert an ld128 number to hex; otherwise, I would include a comment with the hex representation. > The double precision pi*hi and pi*lo also seemed to be sloppily > rounded. My version has less precision (28 low zero bits in pi_hi > instead of the minimal 24), while yours has 27 low zero bits in > pihi. The explanation for this is something like: > - I copied pi_hi from another file with slightly different requirements > - the choice makes pi_lo positive but pilo negative. The other file > might have need pi_lo positive, or just extra zero bits in pi_hi. > Most likely just the latter. > - your 27 is apparently from the calculation 53 - 53 / 2 for almost > even splitting of x, but the splitting of x is uneven (24 + 29). The last one is the closest. I did p/2+1. Yes, the splitting for x is too sloppy. -- Steve 20170425 https://www.youtube.com/watch?v=VWUpyCsUKR4 20161221 https://www.youtube.com/watch?v=IbCHE-hONow