From owner-freebsd-numerics@freebsd.org Fri Mar 8 16:24:37 2019 Return-Path: Delivered-To: freebsd-numerics@mailman.ysv.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1]) by mailman.ysv.freebsd.org (Postfix) with ESMTP id 457EC150F484 for ; Fri, 8 Mar 2019 16:24:37 +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.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits) server-signature RSA-PSS (4096 bits) client-signature RSA-PSS (2048 bits) client-digest SHA256) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 025A6820A8 for ; Fri, 8 Mar 2019 16:24:35 +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 x28GOWLb031509 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO); Fri, 8 Mar 2019 08:24:32 -0800 (PST) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.15.2/8.15.2/Submit) id x28GOWtr031508; Fri, 8 Mar 2019 08:24:32 -0800 (PST) (envelope-from sgk) Date: Fri, 8 Mar 2019 08:24:32 -0800 From: Steve Kargl To: Bruce Evans Cc: freebsd-numerics@freebsd.org Subject: Re: cexpl() (was: Re: Update ENTERI() macro) Message-ID: <20190308162432.GA31392@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20190306055201.GA40298@troutmask.apl.washington.edu> <20190306225811.P2731@besplex.bde.org> <20190306184829.GA44023@troutmask.apl.washington.edu> <20190307061214.R4911@besplex.bde.org> <20190306214233.GA23159@troutmask.apl.washington.edu> <20190307144315.N932@besplex.bde.org> <20190307044447.GA16298@troutmask.apl.washington.edu> <20190307163958.V1293@besplex.bde.org> <20190307195436.GA20856@troutmask.apl.washington.edu> <20190308225807.D6410@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20190308225807.D6410@besplex.bde.org> User-Agent: Mutt/1.11.2 (2019-01-07) X-Rspamd-Queue-Id: 025A6820A8 X-Spamd-Bar: + Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [1.13 / 15.00]; ARC_NA(0.00)[]; HAS_REPLYTO(0.00)[sgk@troutmask.apl.washington.edu]; FROM_HAS_DN(0.00)[]; TO_DN_SOME(0.00)[]; NEURAL_SPAM_SHORT(0.64)[0.636,0]; NEURAL_HAM_LONG(-0.73)[-0.728,0]; MIME_GOOD(-0.10)[text/plain]; RCVD_TLS_LAST(0.00)[]; DMARC_NA(0.00)[washington.edu]; AUTH_NA(1.00)[]; REPLYTO_ADDR_EQ_FROM(0.00)[]; RCVD_COUNT_THREE(0.00)[3]; TO_MATCH_ENVRCPT_SOME(0.00)[]; RCVD_IN_DNSWL_MED(-0.20)[21.76.95.128.list.dnswl.org : 127.0.11.2]; MX_GOOD(-0.01)[cached: troutmask.apl.washington.edu]; RCPT_COUNT_TWO(0.00)[2]; NEURAL_SPAM_MEDIUM(0.49)[0.490,0]; R_SPF_NA(0.00)[]; FREEMAIL_TO(0.00)[optusnet.com.au]; FROM_EQ_ENVFROM(0.00)[]; R_DKIM_NA(0.00)[]; MIME_TRACE(0.00)[0:+]; ASN(0.00)[asn:73, ipnet:128.95.0.0/16, country:US]; MID_RHS_MATCH_FROM(0.00)[]; IP_SCORE(0.05)[ip: (0.10), ipnet: 128.95.0.0/16(0.15), asn: 73(0.05), country: US(-0.07)] X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.29 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: Fri, 08 Mar 2019 16:24:37 -0000 On Sat, Mar 09, 2019 at 12:51:07AM +1100, Bruce Evans wrote: > On Thu, 7 Mar 2019, Steve Kargl wrote: > > > On Thu, Mar 07, 2019 at 05:36:22PM +1100, Bruce Evans wrote: > >> On Wed, 6 Mar 2019, Steve Kargl wrote: > >> > >>> On Thu, Mar 07, 2019 at 03:09:06PM +1100, Bruce Evans wrote: > >* ... > >>> /* > >>> * exp_ovfl = 11356.5234062941439497, 0x400c, 0xb17217f7d1cf79ac > >>> * cexp_ovfl = 22755.3287906024445633, 0x400d, 0xb1c6a8573de9768c > >>> */ > > I checked that these are correct. Any rounding should be down for > exp_ovfl and up for cexp_ovfl. > > >>> if ((hx == 0x400c && lx > 0xb17217f7d1cf79acULL) || > >>> (hx == 0x400d && lx < 0xb1c6a8573de9768cULL)) { > >> > >> Don't duplicate the hex constants in the comment. The limits can be fuzzy > >> so don't need so many decimal digits or nonzero nonzero bits or LD80C() > >> to round them (s_cexp.c uses only 1 decimal digit and doesn't test the > >> lower 32 bits in the mantissa; here we have more bits than we need in lx). > > > > Being fuzzy was the movitation for the new macro, I was only using > > the upper 32 bits of lx in the original cexpl() I posted. > > IIRC, that was done by extracting into a 32-bit lx. This could be written > as (lx >> 32) > 0xb17217f7 (or a fuzzier value), and it might be better > to do that anyway. This is what I do in my WIP ccosh. There are 3 thresholds of ~22.9, ~11356, and ~22755. I'm casting the result to uint32_t, so something like (uint32_t)(lx >> 32) > 0xb17000000. Do I need the cast and/or do you prefer it? > >> Maybe not use bit-fiddling at all for this (in other precisions too). > >> The more important exp() functions check thresholds in FP. Here x can > >> be NaN so FP comparisons with it don't work right. That seems to be > >> the only reason to check in bits here. > >> > >>> /* > >>> * x is between exp_ovfl and cexp_ovfl, so we must scale to > >>> * avoid overflow in exp(x). > >>> */ > >>> RETURNI(__ldexp_cexpl(z, 1)); > >> > >> Did your tests cover this case, and does it use the fixed version of the > >> old __ldexp_cexpl() in k_expl.h? > > > > No. The test I reported was restricted to -11350 < x < 11350. > > Spot checked x outside that range. I had the conditional wrong > > for x < -exp_ovfl. Also, the comment is actaully a little optimistic. > > One cannot avoid overflow. For example, x = 11357 is just above > > exp_ovfl, and exp(x) = inf. __ldexp_cexpl() does scaling and exploits > > |cos(y)| <= 1 and |sin(y)| <= 1 to try to extend the range of cexpl(). > > I understand this better now. __ldexp_cexpl() doesn't work for all large > x, so callers must only use it for a range of values above cexp_ovl. > The correctness of this is unobvious. It depends on sin(y) and cos(y) > never underflowing to 0. Otherwise, this method would give Inf * 0 = NaN > for finite x and y, and the scaling method would give finite * 0 * 2**k = 0. cos(y) can never be 0 and sin(y) is only 0 at y = 0 (see sinpi, cospi implementations). y=0 is a special case and handled separately. We can have Inf*subnormal when y is near 0. -- Steve