From owner-freebsd-numerics@freebsd.org Sat Feb 2 00:25:01 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 31D3014C68AF for ; Sat, 2 Feb 2019 00:25:01 +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 1AC668B581 for ; Sat, 2 Feb 2019 00:25:00 +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 x120OrkF017782 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Fri, 1 Feb 2019 16:24:53 -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 x120OrJu017781 for freebsd-numerics@freebsd.org; Fri, 1 Feb 2019 16:24:53 -0800 (PST) (envelope-from sgk) Date: Fri, 1 Feb 2019 16:24:53 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: cexp[f] patch Message-ID: <20190202002453.GA17741@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.11.2 (2019-01-07) X-Rspamd-Queue-Id: 1AC668B581 X-Spamd-Bar: +++ Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [3.58 / 15.00]; ARC_NA(0.00)[]; HAS_REPLYTO(0.00)[sgk@troutmask.apl.washington.edu]; FROM_HAS_DN(0.00)[]; TO_MATCH_ENVRCPT_ALL(0.00)[]; NEURAL_SPAM_MEDIUM(0.95)[0.953,0]; MIME_GOOD(-0.10)[text/plain]; PREVIOUSLY_DELIVERED(0.00)[freebsd-numerics@freebsd.org]; TO_DN_NONE(0.00)[]; AUTH_NA(1.00)[]; RCPT_COUNT_ONE(0.00)[1]; RCVD_COUNT_THREE(0.00)[3]; RCVD_TLS_LAST(0.00)[]; NEURAL_SPAM_SHORT(0.84)[0.839,0]; RCVD_IN_DNSWL_MED(-0.20)[21.76.95.128.list.dnswl.org : 127.0.11.2]; MX_GOOD(-0.01)[troutmask.apl.washington.edu]; DMARC_NA(0.00)[washington.edu]; NEURAL_SPAM_LONG(0.95)[0.949,0]; R_SPF_NA(0.00)[]; REPLYTO_ADDR_EQ_FROM(0.00)[]; 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.15)[ip: (0.36), ipnet: 128.95.0.0/16(0.39), asn: 73(0.09), 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: Sat, 02 Feb 2019 00:25:01 -0000 The following patch utilizes sincos[f] in the computation of cexp[f]. For 20 million random z=x+Iy drawn in the box defined by x,y in [0,MAX_EXP*LN2], this amounts to a 8.7% and 11.4% speed improvement over computing sin[f] and cos[f] individually. Index: msun/src/s_cexp.c =================================================================== --- msun/src/s_cexp.c (revision 2135) +++ msun/src/s_cexp.c (working copy) @@ -41,7 +41,7 @@ double complex cexp(double complex z) { - double x, y, exp_x; + double c, exp_x, s, x, y; uint32_t hx, hy, lx, ly; x = creal(z); @@ -55,8 +55,10 @@ return (CMPLX(exp(x), y)); EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if (((hx & 0x7fffffff) | lx) == 0) - return (CMPLX(cos(y), sin(y))); + if (((hx & 0x7fffffff) | lx) == 0) { + sincos(y, &s, &c); + return (CMPLX(c, s)); + } if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { @@ -86,6 +88,7 @@ * - x = NaN (spurious inexact exception from y) */ exp_x = exp(x); - return (CMPLX(exp_x * cos(y), exp_x * sin(y))); + sincos(y, &s, &c); + return (CMPLX(exp_x * c, exp_x * s)); } } Index: msun/src/s_cexpf.c =================================================================== --- msun/src/s_cexpf.c (revision 2135) +++ msun/src/s_cexpf.c (working copy) @@ -41,7 +41,7 @@ float complex cexpf(float complex z) { - float x, y, exp_x; + float c, exp_x, s, x, y; uint32_t hx, hy; x = crealf(z); @@ -55,8 +55,10 @@ return (CMPLXF(expf(x), y)); GET_FLOAT_WORD(hx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if ((hx & 0x7fffffff) == 0) - return (CMPLXF(cosf(y), sinf(y))); + if ((hx & 0x7fffffff) == 0) { + sincosf(y, &s, &c); + return (CMPLXF(c, s)); + } if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { @@ -86,6 +88,7 @@ * - x = NaN (spurious inexact exception from y) */ exp_x = expf(x); - return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y))); + sincosf(y, &s, &c); + return (CMPLXF(exp_x * c, exp_x * s)); } } -- Steve