From owner-freebsd-numerics@freebsd.org Thu Mar 7 04:44:53 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 A1F5E150FABC for ; Thu, 7 Mar 2019 04:44:53 +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 13C418D285 for ; Thu, 7 Mar 2019 04:44:51 +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 x274ilAZ016396 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO); Wed, 6 Mar 2019 20:44:47 -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 x274ilQo016395; Wed, 6 Mar 2019 20:44:47 -0800 (PST) (envelope-from sgk) Date: Wed, 6 Mar 2019 20:44:47 -0800 From: Steve Kargl To: Bruce Evans Cc: freebsd-numerics@freebsd.org Subject: Re: Update ENTERI() macro Message-ID: <20190307044447.GA16298@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20190227161906.GA77785@troutmask.apl.washington.edu> <20190228060920.R4413@besplex.bde.org> <20190304212159.GA12587@troutmask.apl.washington.edu> <20190305153243.Y1349@besplex.bde.org> <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> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20190307144315.N932@besplex.bde.org> User-Agent: Mutt/1.11.2 (2019-01-07) X-Rspamd-Queue-Id: 13C418D285 X-Spamd-Bar: + Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [1.05 / 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.96)[0.957,0]; NEURAL_HAM_LONG(-0.72)[-0.719,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.08)[0.079,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: Thu, 07 Mar 2019 04:44:54 -0000 On Thu, Mar 07, 2019 at 03:09:06PM +1100, Bruce Evans wrote: > On Wed, 6 Mar 2019, Steve Kargl wrote: > > > On Thu, Mar 07, 2019 at 06:30:42AM +1100, Bruce Evans wrote: > >> ... > >> I now see that you implemented 2 more versions of __ldexpl_cexpl() by > >> cloning the old double precision version. Apparently the includes > >> are to unpolluted for the compiler to see the multiple versions :-). > >> > >> Using the version in k_expl.h almost forces inlining of expl()'s kernel > >> and its large tables, just like for hyperbolic functions. This wastes > >> a lot of space, especially for duplicating the tables. It is only a > >> small optimization for time. It is done for the hyperbolic functions > >> to get this optimization, and for __ldexpl_cexpl() just for convenience. > > > > The version in k_expl.h has 2 bugs. You note the first (cos instead > > of cosl). The second is > > > > In file included from /data/kargl/trunk/math/libm/msun/ld80/s_cexpl.c:43: > > /data/kargl/trunk/math/libm/msun/ld80/k_expl.h:288:22: error: magnitude of > > floating-point constant too large for type 'double'; maximum is > > 1.7976931348623157E+308 [-Werror,-Wliteral-range] > > exp_x = (lo + hi) * 0x1p16382; > > ^ > > It is missing an L suffix. Clearly not tested. > > >>> [errors for cexpl()] > >> > >> Check a few denormals, Infs and NaNs. > > > > Exceptional cases give > > > > % ./testl -e > > cexpl(1, inf) = (nan,nan) expecting (nan,nan) > > cexpl(1,-inf) = (nan,nan) expecting (nan,nan) > > cexpl(nan,-inf) = (nan,nan) expecting (nan,nan) > > cexpl(nan, inf) = (nan,nan) expecting (nan,nan) > > cexpl(inf, inf) = (inf,nan) expecting (inf,nan) > > cexpl(inf,-inf) = (inf,nan) expecting (inf,nan) > > cexpl(inf,nan) = (inf,nan) expecting (inf,nan) > > cexpl(-inf,nan) = (0,0) expecting (0,0) > > cexpl(-inf,-inf) = (0,0) expecting (0,0) > > cexpl(-inf, inf) = (0,0) expecting (0,0) > > But does it preserve NaN bits as much as possible, and produce the same > NaN bits as double precision after conversion to double precision? The > report spells NaN and Inf using C99's fuzzy bad names. Even the C99 > specification spells NaN and Inf correctly in Appendix F. > > The complex hyperbolic functions have better comments, so are almost > an easier place to start. > I suppose it preserves the NaN bits as well as cexpf and cexp as it uses the same algorthim. /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * * Copyright (c) 2011 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. * * src/s_cexp.c converted to long double complex by Steven G. Kargl */ #include __FBSDID("$FreeBSD$"); #include #include #ifdef __i386__ #include #endif #include "fpmath.h" #include "math.h" #include "math_private.h" #include "k_expl.h" long double complex cexpl (long double complex z) { long double c, exp_x, s, x, y; uint64_t lx, ly; uint16_t hx, hy, ix, iy; ENTERI(); x = creall(z); y = cimagl(z); EXTRACT_LDBL80_WORDS(hx, lx, x); EXTRACT_LDBL80_WORDS(hy, ly, y); ix = hx&0x7fff; iy = hy&0x7fff; /* cexp(0 + I y) = cos(y) + I sin(y) */ if ((ix | lx) == 0) { sincosl(y, &s, &c); RETURNI(CMPLXL(c, s)); } /* cexp(x + I 0) = exp(x) + I 0 */ if ((iy | ly) == 0) RETURNI(CMPLXL(expl(x), y)); if (iy == 0x7fff) { if (ix < 0x7fff || (ix == 0x7fff && (lx & 0x7fffffffffffffffULL) != 0)) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ RETURNI(CMPLXL(y - y, y - y)); } else if (hx & 0x8000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ RETURNI(CMPLXL(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ RETURNI(CMPLXL(x, y - y)); } } /* * exp_ovfl = 11356.5234062941439497, 0x400c, 0xb17217f7d1cf79ac * cexp_ovfl = 22755.3287906024445633, 0x400d, 0xb1c6a8573de9768c */ if ((hx == 0x400c && lx > 0xb17217f7d1cf79acULL) || (hx == 0x400d && lx < 0xb1c6a8573de9768cULL)) { /* * x is between exp_ovfl and cexp_ovfl, so we must scale to * avoid overflow in exp(x). */ RETURNI(__ldexp_cexpl(z, 1)); } else { /* * Cases covered here: * - x < exp_ovfl and exp(x) won't overflow (common case) * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 * - x = +-Inf (generated by exp()) * - x = NaN (spurious inexact exception from y) */ exp_x = expl(x); sincosl(y, &s, &c); RETURNI(CMPLXL(exp_x * c, exp_x * s)); } } -- Steve 20170425 https://www.youtube.com/watch?v=VWUpyCsUKR4 20161221 https://www.youtube.com/watch?v=IbCHE-hONow