From owner-freebsd-numerics@freebsd.org Sun Feb 3 17:31:06 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 3460514AC74A for ; Sun, 3 Feb 2019 17:31:06 +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 7DDB5778EA for ; Sun, 3 Feb 2019 17:31:04 +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 x13HUuY2042514 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Sun, 3 Feb 2019 09:30:56 -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 x13HUu8f042513 for freebsd-numerics@freebsd.org; Sun, 3 Feb 2019 09:30:56 -0800 (PST) (envelope-from sgk) Date: Sun, 3 Feb 2019 09:30:56 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Implement cexpl Message-ID: <20190203173056.GA42499@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: 7DDB5778EA X-Spamd-Bar: +++ Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [3.57 / 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.91)[0.906,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.94)[0.941,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.88)[0.884,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.35), ipnet: 128.95.0.0/16(0.38), 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: Sun, 03 Feb 2019 17:31:06 -0000 The following patch implements cexpl() for ld80 and ld128 architectures. The kernel files for large x are named ld{80,128}/k_cexpl.c to avoid confusion with the ld{80,128}/k_expl.h files. I do not have access to ld128, so the ld128 does not use bit twiddling. I cannot test the ld128 implementation, but I believe that it works (for some definition of work). The large patch does not contain the needed change to include/complex.h. Index: include/complex.h =================================================================== --- include/complex.h (revision 343477) +++ include/complex.h (working copy) @@ -98,6 +98,8 @@ double complex ccosh(double complex); float complex ccoshf(float complex); double complex cexp(double complex); float complex cexpf(float complex); +long double complex + cexpl(long double complex); double cimag(double complex) __pure2; float cimagf(float complex) __pure2; long double cimagl(long double complex) __pure2; Here's a svn logfile. * Makefile: . Add k_cexpl.c and s_cexpl.c to build. * src/math_private.h: . Add ENTERC and RETURNC to allow toggling of FP_PE in i686-class hardware for long double complex. . Add prototype for __ldexp_cexpl(). * src/s_cexp.c: . Include float.h to get LDBL_MANT_DIG so that a weak reference can be add for LDBL_MANT_DIG == 53 targets. . Flip the y == 0 and x == 0 cases. This is in anticipation of allowing GCC to transform CMPLX(cos(x), sin(x)) into cexp(CMPLX(0,x)), which can then be mapped to the sincos() builtin. * src/s_cexpf.c: . Flip the y == 0 and x == 0 cases. This is in anticipation of allowing GCC to transform CMPLX(cos(x), sin(x)) into cexp(CMPLX(0,x)), which can then be mapped to the sincos() builtin. * src/k_exp.c: . Declare c and s, and re-order declarations. . Use sincos(). * src/k_expf.c: . Declare c and s, and re-order declarations. . Use sincosf(). * ld80/k_cexpl.c: . Implementation of long double complex kernels for cexpl() for x large. . Conversion of src/k_cexp.c from double complex to long double complex. * ld80/s_cexpl.c: . Implementation of long double complex cexpl(). . Conversion of src/k_cexp.c from double complex to long double complex. * ld128/k_cexpl.c: . FIXME: I do not have access to ld128 hardware. Someone who cares about ld128 needs to fix this implementation of long double complex kernels for cexpl() for x large. . The functions defined here are currently unused. * ld128/s_cexpl.c: . FIXME: I do not have access to ld128 hardware. Someone who cares about ld128 needs to fix this implementation of long double complex cexpl() to use bit twiddling. It currently uses isinf() and isnan() where needed and floating point comparisons. And now the diff. Index: Makefile =================================================================== --- Makefile (revision 2141) +++ Makefile (working copy) @@ -133,7 +133,9 @@ COMMON_SRCS+= catrigl.c \ e_lgammal.c e_lgammal_r.c e_powl.c \ e_remainderl.c e_sinhl.c e_sqrtl.c \ invtrig.c k_cosl.c k_sinl.c k_tanl.c \ + k_cexpl.c \ s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c \ + s_cexpl.c \ s_clogl.c s_cosl.c s_cospil.c s_cprojl.c \ s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \ s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \ Index: src/math_private.h =================================================================== --- src/math_private.h (revision 2141) +++ src/math_private.h (working copy) @@ -357,9 +357,24 @@ do { \ fpsetprec(__oprec); \ return; \ } while (0) +#define ENTERC() ENTERCT(long double complex) +#define ENTERCT(returntype) \ + returntype __retval; \ + fp_prec_t __oprec; \ + \ + if ((__oprec = fpgetprec()) != FP_PE) \ + fpsetprec(FP_PE) +#define RETURNC(x) do { \ + __retval = (x); \ + if (__oprec != FP_PE) \ + fpsetprec(__oprec); \ + RETURNF(__retval); \ +} while (0) #else +#define ENTERC(x) #define ENTERI(x) #define ENTERIT(x) +#define RETURNC(x) RETURNF(x) #define RETURNI(x) RETURNF(x) #define ENTERV() #define RETURNV() return @@ -919,6 +934,10 @@ float complex __ldexp_cexpf(float complex,int); long double __kernel_sinl(long double, long double, int); long double __kernel_cosl(long double, long double); long double __kernel_tanl(long double, long double, int); +long double __ldexp_expl(long double, int); +#ifdef _COMPLEX_H +long double complex __ldexp_cexpl(long double complex, int); +#endif #include "math_sgk.h" Index: src/s_cexp.c =================================================================== --- src/s_cexp.c (revision 2141) +++ src/s_cexp.c (working copy) @@ -30,6 +30,7 @@ __FBSDID("$FreeBSD: head/lib/msun/src/s_cexp.c 326219 2017-11-26 02:00:33Z pfg $"); #include +#include #include #include "math_private.h" @@ -47,12 +48,6 @@ cexp(double complex z) x = creal(z); y = cimag(z); - EXTRACT_WORDS(hy, ly, y); - hy &= 0x7fffffff; - - /* cexp(x + I 0) = exp(x) + I 0 */ - if ((hy | ly) == 0) - return (CMPLX(exp(x), y)); EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if (((hx & 0x7fffffff) | lx) == 0) { @@ -60,6 +55,13 @@ cexp(double complex z) return (CMPLX(c, s)); } + EXTRACT_WORDS(hy, ly, y); + hy &= 0x7fffffff; + + /* cexp(x + I 0) = exp(x) + I 0 */ + if ((hy | ly) == 0) + return (CMPLX(exp(x), y)); + if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ @@ -92,3 +94,7 @@ cexp(double complex z) return (CMPLX(exp_x * c, exp_x * s)); } } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(cexp, cexpl); +#endif Index: src/s_cexpf.c =================================================================== --- src/s_cexpf.c (revision 2141) +++ src/s_cexpf.c (working copy) @@ -47,18 +47,19 @@ cexpf(float complex z) x = crealf(z); y = cimagf(z); - GET_FLOAT_WORD(hy, y); - hy &= 0x7fffffff; - - /* cexp(x + I 0) = exp(x) + I 0 */ - if (hy == 0) - return (CMPLXF(expf(x), y)); GET_FLOAT_WORD(hx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if ((hx & 0x7fffffff) == 0) { sincosf(y, &s, &c); return (CMPLXF(c, s)); } + + GET_FLOAT_WORD(hy, y); + hy &= 0x7fffffff; + + /* cexp(x + I 0) = exp(x) + I 0 */ + if (hy == 0) + return (CMPLXF(expf(x), y)); if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { Index: src/k_exp.c =================================================================== --- src/k_exp.c (revision 2141) +++ src/k_exp.c (working copy) @@ -88,7 +88,7 @@ __ldexp_exp(double x, int expt) double complex __ldexp_cexp(double complex z, int expt) { - double x, y, exp_x, scale1, scale2; + double c, exp_x, s, scale1, scale2, x, y; int ex_expt, half_expt; x = creal(z); @@ -105,6 +105,7 @@ __ldexp_cexp(double complex z, int expt) half_expt = expt - half_expt; INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); - return (CMPLX(cos(y) * exp_x * scale1 * scale2, - sin(y) * exp_x * scale1 * scale2)); + sincos(y, &s, &c); + return (CMPLX(c * exp_x * scale1 * scale2, + s * exp_x * scale1 * scale2)); } Index: src/k_expf.c =================================================================== --- src/k_expf.c (revision 2141) +++ src/k_expf.c (working copy) @@ -71,7 +71,7 @@ __ldexp_expf(float x, int expt) float complex __ldexp_cexpf(float complex z, int expt) { - float x, y, exp_x, scale1, scale2; + float c, exp_x, s, scale1, scale2, x, y; int ex_expt, half_expt; x = crealf(z); @@ -84,6 +84,7 @@ __ldexp_cexpf(float complex z, int expt) half_expt = expt - half_expt; SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); - return (CMPLXF(cosf(y) * exp_x * scale1 * scale2, - sinf(y) * exp_x * scale1 * scale2)); + sincosf(y, &s, &c); + return (CMPLXF(c * exp_x * scale1 * scale2, + s * exp_x * scale1 * scale2)); } Index: ld80/k_cexpl.c =================================================================== --- ld80/k_cexpl.c (revision 2141) +++ ld80/k_cexpl.c (working copy) @@ -24,31 +24,34 @@ * 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/k_exp.c converted to long double complex by Steven G. Kargl */ #include -__FBSDID("$FreeBSD: head/lib/msun/src/k_exp.c 326219 2017-11-26 02:00:33Z pfg $"); +__FBSDID("$FreeBSD$"); #include +#include "fpmath.h" #include "math.h" #include "math_private.h" -static const uint32_t k = 1799; /* constant for reduction */ -static const double kln2 = 1246.97177782734161156; /* k * ln2 */ - +static const uint32_t k = 22956; /* constant for reduction */ +static const union IEEEl2bits +kln2u = LD80C(0xf89f8bf509c862ca, 13, 1.59118866769341045231e+04L); +#define kln2 (kln2u.e) /* * Compute exp(x), scaled to avoid spurious overflow. An exponent is * returned separately in 'expt'. * - * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 - * Output: 2**1023 <= y < 2**1024 + * Input: ln(LDBL_MAX) <= x < ln(2 * LDBL_MAX / LDBL_MIN_DENORM) ~= 22756.0 + * Output: 2**16383 <= y < 2**16384 */ -static double -__frexp_exp(double x, int *expt) +static long double +__frexp_expl(long double x, int *expt) { - double exp_x; - uint32_t hx; + union IEEEl2bits exp_x; /* * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to @@ -56,11 +59,10 @@ __frexp_exp(double x, int *expt) * exp_x to MAX_EXP so that the result can be multiplied by * a tiny number without losing accuracy due to denormalization. */ - exp_x = exp(x - kln2); - GET_HIGH_WORD(hx, exp_x); - *expt = (hx >> 20) - (0x3ff + 1023) + k; - SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); - return (exp_x); + exp_x.e = expl(x - kln2); + *expt = exp_x.bits.exp - (0x3fff + 16383) + k - 1; + exp_x.bits.exp = 0x3fff + 16383; + return (exp_x.e); } /* @@ -73,27 +75,30 @@ __frexp_exp(double x, int *expt) * has filtered out very large x, for which overflow would be inevitable. */ -double -__ldexp_exp(double x, int expt) +long double +__ldexp_expl(long double x, int expt) { - double exp_x, scale; + union IEEEl2bits scale; + long double exp_x; int ex_expt; - exp_x = __frexp_exp(x, &ex_expt); + exp_x = __frexp_expl(x, &ex_expt); expt += ex_expt; - INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); - return (exp_x * scale); + scale.e = 1; + scale.bits.exp = 0x3fff + expt; + return (exp_x * scale.e); } -double complex -__ldexp_cexp(double complex z, int expt) +long double complex +__ldexp_cexpl(long double complex z, int expt) { - double x, y, exp_x, scale1, scale2; + union IEEEl2bits scale1, scale2; + long double c, exp_x, s, x, y; int ex_expt, half_expt; - x = creal(z); - y = cimag(z); - exp_x = __frexp_exp(x, &ex_expt); + x = creall(z); + y = cimagl(z); + exp_x = __frexp_expl(x, &ex_expt); expt += ex_expt; /* @@ -101,10 +106,13 @@ __ldexp_cexp(double complex z, int expt) * compensate for scalbn being horrendously slow. */ half_expt = expt / 2; - INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); + scale1.e = 1; + scale1.bits.exp = 0x3fff + half_expt; half_expt = expt - half_expt; - INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); + scale2.e = 1; + scale2.bits.exp = 0x3fff + half_expt + 1; - return (CMPLX(cos(y) * exp_x * scale1 * scale2, - sin(y) * exp_x * scale1 * scale2)); + sincosl(y, &s, &c); + return (CMPLXL(c * exp_x * scale1.e * scale2.e, + s * exp_x * scale1.e * scale2.e)); } Index: ld80/s_cexpl.c =================================================================== --- ld80/s_cexpl.c (revision 2141) +++ ld80/s_cexpl.c (working copy) @@ -24,61 +24,72 @@ * 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: head/lib/msun/src/s_cexp.c 326219 2017-11-26 02:00:33Z pfg $"); +__FBSDID("$FreeBSD$"); #include +#include #include +#ifdef __i386__ +#include +#endif + +#include "fpmath.h" #include "math_private.h" static const uint32_t -exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ -cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ +exp_ovfl = 0xb17217f7, /* high bits of MAX_EXP * ln2 ~= 11356 */ +cexp_ovfl = 0xb1c6a857; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ -double complex -cexp(double complex z) +long double complex +cexpl (long double complex z) { - double c, exp_x, s, x, y; - uint32_t hx, hy, lx, ly; + union IEEEl2bits ux, uy; + long double c, exp_x, s, x, y; + uint16_t hx, hy; - x = creal(z); - y = cimag(z); + ENTERC(); - EXTRACT_WORDS(hy, ly, y); - hy &= 0x7fffffff; + x = ux.e = creall(z); + y = uy.e = cimagl(z); - /* cexp(x + I 0) = exp(x) + I 0 */ - if ((hy | ly) == 0) - return (CMPLX(exp(x), y)); - EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if (((hx & 0x7fffffff) | lx) == 0) { - sincos(y, &s, &c); - return (CMPLX(c, s)); + hx = ux.xbits.expsign; + if ((ux.bits.manh | ux.bits.manl | (hx & 0x7fff)) == 0) { + sincosl(y, &s, &c); + RETURNC(CMPLXL(c, s)); } - if (hy >= 0x7ff00000) { - if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { + /* cexp(x + I 0) = exp(x) + I 0 */ + hy = uy.xbits.expsign; + if ((uy.bits.manh | uy.bits.manl | (hy & 0x7fff)) == 0) + RETURNC(CMPLXL(expl(x), y)); + + if (hy >= 0x7fff) { + if (hx < 0x7fff || + (hx == 0x7fff && (ux.bits.manh & 0x7fffffff) != 0)) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (CMPLX(y - y, y - y)); - } else if (hx & 0x80000000) { + RETURNC(CMPLXL(y - y, y - y)); + } else if (hx & 0x8000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (CMPLX(0.0, 0.0)); + RETURNC(CMPLXL(0.0L, 0.0L)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (CMPLX(x, y - y)); + RETURNC(CMPLXL(x, y - y)); } } - if (hx >= exp_ovfl && hx <= cexp_ovfl) { + if (hx >= 0x400c && hx <= 0x400d) { /* - * x is between 709.7 and 1454.3, so we must scale to avoid + * x is between 11356 and 22755, so we must scale to avoid * overflow in exp(x). */ - return (__ldexp_cexp(z, 0)); + RETURNC(__ldexp_cexpl(z, 0)); } else { /* * Cases covered here: @@ -87,8 +98,8 @@ cexp(double complex z) * - x = +-Inf (generated by exp()) * - x = NaN (spurious inexact exception from y) */ - exp_x = exp(x); - sincos(y, &s, &c); - return (CMPLX(exp_x * c, exp_x * s)); + exp_x = expl(x); + sincosl(y, &s, &c); + RETURNC(CMPLXL(exp_x * c, exp_x * s)); } } Index: ld128/k_cexpl.c =================================================================== --- ld128/k_cexpl.c (revision 2141) +++ ld128/k_cexpl.c (working copy) @@ -34,8 +34,14 @@ __FBSDID("$FreeBSD: head/lib/msun/src/k_exp.c 326219 2 #include "math.h" #include "math_private.h" +#warning libm is broken on ld128 architectures +#warning functions defined here are currently unused +#warning someone who cares needs to convert src/k_exp.c + +#if 0 static const uint32_t k = 1799; /* constant for reduction */ static const double kln2 = 1246.97177782734161156; /* k * ln2 */ +#endif /* * Compute exp(x), scaled to avoid spurious overflow. An exponent is @@ -44,9 +50,10 @@ static const double kln2 = 1246.97177782734161156; /* * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 * Output: 2**1023 <= y < 2**1024 */ -static double -__frexp_exp(double x, int *expt) +static long double +__frexp_expl(long double x, int *expt) { +#if 0 double exp_x; uint32_t hx; @@ -61,6 +68,8 @@ __frexp_exp(double x, int *expt) *expt = (hx >> 20) - (0x3ff + 1023) + k; SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); return (exp_x); +#endif + return x * expt; } /* @@ -73,9 +82,10 @@ __frexp_exp(double x, int *expt) * has filtered out very large x, for which overflow would be inevitable. */ -double -__ldexp_exp(double x, int expt) +long double +__ldexp_expl(long double x, int expt) { +#if 0 double exp_x, scale; int ex_expt; @@ -83,11 +93,14 @@ __ldexp_exp(double x, int expt) expt += ex_expt; INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); return (exp_x * scale); +#endif + return x * expt; } -double complex -__ldexp_cexp(double complex z, int expt) +long double complex +__ldexp_cexpl(long double complex z, int expt) { +#if 0 double x, y, exp_x, scale1, scale2; int ex_expt, half_expt; @@ -105,6 +118,9 @@ __ldexp_cexp(double complex z, int expt) half_expt = expt - half_expt; INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); - return (CMPLX(cos(y) * exp_x * scale1 * scale2, - sin(y) * exp_x * scale1 * scale2)); + sincos(y, &s, &c); + return (CMPLX(c) * exp_x * scale1 * scale2, + s * exp_x * scale1 * scale2)); +#endif + return z * expt; } Index: ld128/s_cexpl.c =================================================================== --- ld128/s_cexpl.c (revision 2141) +++ ld128/s_cexpl.c (working copy) @@ -1,7 +1,7 @@ /*- * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * - * Copyright (c) 2011 David Schultz + * Copyright (c) 2019 Steven G. Kargl * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -27,68 +27,49 @@ */ #include -__FBSDID("$FreeBSD: head/lib/msun/src/s_cexp.c 326219 2017-11-26 02:00:33Z pfg $"); +__FBSDID("$FreeBSD$"); #include +#include #include #include "math_private.h" -static const uint32_t -exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ -cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ +#warning libm is broken on ld128 architectures +#warning someone who cares needs to convert src/s_cexp.c -double complex -cexp(double complex z) +long double complex +cexpl(long double complex z) { - double c, exp_x, s, x, y; - uint32_t hx, hy, lx, ly; + long double c, exp_x, s, x, y; - x = creal(z); - y = cimag(z); + x = creall(z); + y = cimagl(z); - EXTRACT_WORDS(hy, ly, y); - hy &= 0x7fffffff; - - /* cexp(x + I 0) = exp(x) + I 0 */ - if ((hy | ly) == 0) - return (CMPLX(exp(x), y)); - EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ - if (((hx & 0x7fffffff) | lx) == 0) { - sincos(y, &s, &c); - return (CMPLX(c, s)); + if (x == 0) { + sincosl(y, &s, &c); + return (CMPLXL(c, s)); } - if (hy >= 0x7ff00000) { - if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { + /* cexp(x + I 0) = exp(x) + I 0 */ + if (y == 0) + return (CMPLXL(expl(x), y)); + + if (!isfinite(y)) { + if (isfinite(x) || isnan(x)) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (CMPLX(y - y, y - y)); - } else if (hx & 0x80000000) { + return (CMPLXL(y - y, y - y)); + } else if (isinf(x) && copysignl(1.L, x) < 0) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (CMPLX(0.0, 0.0)); + return (CMPLXL(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (CMPLX(x, y - y)); + return (CMPLXL(x, y - y)); } } - if (hx >= exp_ovfl && hx <= cexp_ovfl) { - /* - * x is between 709.7 and 1454.3, so we must scale to avoid - * overflow in exp(x). - */ - return (__ldexp_cexp(z, 0)); - } 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 = exp(x); - sincos(y, &s, &c); - return (CMPLX(exp_x * c, exp_x * s)); - } + exp_x = expl(x); + sincosl(y, &s, &c); + return (CMPLXL(exp_x * c, exp_x * s)); } -- Steve From owner-freebsd-numerics@freebsd.org Sun Feb 3 17:43:03 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 BAD3014ACC63 for ; Sun, 3 Feb 2019 17:43:03 +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 81B2877DAC for ; Sun, 3 Feb 2019 17:43:01 +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 x13Hgxb9042575 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Sun, 3 Feb 2019 09:42:59 -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 x13Hgxtw042574 for freebsd-numerics@freebsd.org; Sun, 3 Feb 2019 09:42:59 -0800 (PST) (envelope-from sgk) Date: Sun, 3 Feb 2019 09:42:59 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: Implement cexpl Message-ID: <20190203174259.GA42562@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20190203173056.GA42499@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20190203173056.GA42499@troutmask.apl.washington.edu> User-Agent: Mutt/1.11.2 (2019-01-07) X-Rspamd-Queue-Id: 81B2877DAC X-Spamd-Bar: +++ Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [3.59 / 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.90)[0.900,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.95)[0.953,0]; 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]; DMARC_NA(0.00)[washington.edu]; NEURAL_SPAM_LONG(0.90)[0.905,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.34), ipnet: 128.95.0.0/16(0.37), 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: Sun, 03 Feb 2019 17:43:03 -0000 On Sun, Feb 03, 2019 at 09:30:56AM -0800, Steve Kargl wrote: > The following patch implements cexpl() for ld80 and ld128 architectures. > > The kernel files for large x are named ld{80,128}/k_cexpl.c to avoid > confusion with the ld{80,128}/k_expl.h files. > > I do not have access to ld128, so the ld128 does not use bit twiddling. > I cannot test the ld128 implementation, but I believe that it works > (for some definition of work). > I forgot to mention that the initial files under ld{80,128} where created via svn cp src/k_exp.c ld80/k_cexpl.c svn cp src/s_cexp.c ld80/s_cexpl.c svn cp src/k_exp.c ld128/k_cexpl.c svn cp src/s_cexp.c ld128/s_cexpl.c so the diff is not against an empty file; but against the copied files. -- Steve From owner-freebsd-numerics@freebsd.org Sun Feb 3 19:34:50 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 3B81C14B01CB for ; Sun, 3 Feb 2019 19:34:50 +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 E86448338C for ; Sun, 3 Feb 2019 19:34:48 +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 x13JYkKA043330 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Sun, 3 Feb 2019 11:34:46 -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 x13JYkAs043329 for freebsd-numerics@freebsd.org; Sun, 3 Feb 2019 11:34:46 -0800 (PST) (envelope-from sgk) Date: Sun, 3 Feb 2019 11:34:46 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: optimization for Bessel function routines Message-ID: <20190203193446.GA43322@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: E86448338C X-Spamd-Bar: +++ Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [3.53 / 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.88)[0.885,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.93)[0.933,0]; 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]; DMARC_NA(0.00)[washington.edu]; NEURAL_SPAM_LONG(0.87)[0.875,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.14)[ip: (0.33), ipnet: 128.95.0.0/16(0.37), 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: Sun, 03 Feb 2019 19:34:50 -0000 * src/e_j0.c: . Replace calls to sin(x) and cos(x) with single call to sincos(). * src/e_j0f.c: . Replace calls to sinf(x) and cosf(x) with single call to sincosf(). * src/e_j1.c: . Replace calls to sin(x) and cos(x) with single call to sincos(). * src/e_j1f.c: . Replace calls to sinf(x) and cosf(x) with single call to sincosf(). * src/e_jn.c: . Replace calls to sin(x) and cos(x) with single call to sincos(). Index: src/e_j0.c =================================================================== --- src/e_j0.c (revision 343477) +++ src/e_j0.c (working copy) @@ -93,8 +93,7 @@ __ieee754_j0(double x) if(ix>=0x7ff00000) return one/(x*x); x = fabs(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sin(x); - c = cos(x); + sincos(x, &s, &c); ss = s-c; cc = s+c; if(ix<0x7fe00000) { /* Make sure x+x does not overflow. */ @@ -173,8 +172,7 @@ __ieee754_y0(double x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - s = sin(x); - c = cos(x); + sincos(x, &s, &c); ss = s-c; cc = s+c; /* Index: src/e_j0f.c =================================================================== --- src/e_j0f.c (revision 343477) +++ src/e_j0f.c (working copy) @@ -55,8 +55,7 @@ __ieee754_j0f(float x) if(ix>=0x7f800000) return one/(x*x); x = fabsf(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sinf(x); - c = cosf(x); + sincosf(x, &s, &c); ss = s-c; cc = s+c; if(ix<0x7f000000) { /* Make sure x+x does not overflow. */ @@ -128,8 +127,7 @@ __ieee754_y0f(float x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - s = sinf(x); - c = cosf(x); + sincosf(x, &s, &c); ss = s-c; cc = s+c; /* Index: src/e_j1.c =================================================================== --- src/e_j1.c (revision 343477) +++ src/e_j1.c (working copy) @@ -94,8 +94,7 @@ __ieee754_j1(double x) if(ix>=0x7ff00000) return one/x; y = fabs(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sin(y); - c = cos(y); + sincos(y, &s, &c); ss = -s-c; cc = s-c; if(ix<0x7fe00000) { /* make sure y+y not overflow */ @@ -159,8 +158,7 @@ __ieee754_y1(double x) /* y1(x<0) = NaN and raise invalid exception. */ if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sin(x); - c = cos(x); + sincos(x, &s, &c); ss = -s-c; cc = s-c; if(ix<0x7fe00000) { /* make sure x+x not overflow */ Index: src/e_j1f.c =================================================================== --- src/e_j1f.c (revision 343477) +++ src/e_j1f.c (working copy) @@ -56,8 +56,7 @@ __ieee754_j1f(float x) if(ix>=0x7f800000) return one/x; y = fabsf(x); if(ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sinf(y); - c = cosf(y); + sincosf(y, &s, &c); ss = -s-c; cc = s-c; if(ix<0x7f000000) { /* make sure y+y not overflow */ @@ -114,8 +113,7 @@ __ieee754_y1f(float x) if(ix==0) return -one/vzero; if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sinf(x); - c = cosf(x); + sincosf(x, &s, &c); ss = -s-c; cc = s-c; if(ix<0x7f000000) { /* make sure x+x not overflow */ Index: src/e_jn.c =================================================================== --- src/e_jn.c (revision 343477) +++ src/e_jn.c (working copy) @@ -54,7 +54,7 @@ double __ieee754_jn(int n, double x) { int32_t i,hx,ix,lx, sgn; - double a, b, temp, di; + double a, b, c, s, temp, di; double z, w; /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) @@ -91,11 +91,12 @@ __ieee754_jn(int n, double x) * 2 -s+c -c-s * 3 s+c c-s */ + sincos(x, &s, &c); switch(n&3) { - case 0: temp = cos(x)+sin(x); break; - case 1: temp = -cos(x)+sin(x); break; - case 2: temp = -cos(x)-sin(x); break; - case 3: temp = cos(x)-sin(x); break; + case 0: temp = c+s; break; + case 1: temp = -c+s; break; + case 2: temp = -c-s; break; + case 3: temp = c-s; break; } b = invsqrtpi*temp/sqrt(x); } else { @@ -216,7 +217,7 @@ __ieee754_yn(int n, double x) { int32_t i,hx,ix,lx; int32_t sign; - double a, b, temp; + double a, b, c, s, temp; EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; @@ -248,11 +249,12 @@ __ieee754_yn(int n, double x) * 2 -s+c -c-s * 3 s+c c-s */ + sincos(x, &s, &c); switch(n&3) { - case 0: temp = sin(x)-cos(x); break; - case 1: temp = -sin(x)-cos(x); break; - case 2: temp = -sin(x)+cos(x); break; - case 3: temp = sin(x)+cos(x); break; + case 0: temp = s-c; break; + case 1: temp = -s-c; break; + case 2: temp = -s+c; break; + case 3: temp = s+c; break; } b = invsqrtpi*temp/sqrt(x); } else { -- Steve From owner-freebsd-numerics@freebsd.org Mon Feb 4 19:32:51 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 9792614AFD91 for ; Mon, 4 Feb 2019 19:32:51 +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 125946EBD2 for ; Mon, 4 Feb 2019 19:32:50 +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 x14JWlSJ050273 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Mon, 4 Feb 2019 11:32: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 x14JWlrx050272 for freebsd-numerics@freebsd.org; Mon, 4 Feb 2019 11:32:47 -0800 (PST) (envelope-from sgk) Date: Mon, 4 Feb 2019 11:32:47 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: Implement cexpl Message-ID: <20190204193247.GA50265@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20190203173056.GA42499@troutmask.apl.washington.edu> <20190203174259.GA42562@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20190203174259.GA42562@troutmask.apl.washington.edu> User-Agent: Mutt/1.11.2 (2019-01-07) X-Rspamd-Queue-Id: 125946EBD2 X-Spamd-Bar: +++ Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [3.38 / 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.94)[0.939,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.75)[0.745,0]; 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]; DMARC_NA(0.00)[washington.edu]; NEURAL_SPAM_LONG(0.88)[0.879,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.13)[ip: (0.28), ipnet: 128.95.0.0/16(0.34), 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: Mon, 04 Feb 2019 19:32:51 -0000 On Sun, Feb 03, 2019 at 09:42:59AM -0800, Steve Kargl wrote: > On Sun, Feb 03, 2019 at 09:30:56AM -0800, Steve Kargl wrote: > > The following patch implements cexpl() for ld80 and ld128 architectures. > > > > The kernel files for large x are named ld{80,128}/k_cexpl.c to avoid > > confusion with the ld{80,128}/k_expl.h files. > > > > I do not have access to ld128, so the ld128 does not use bit twiddling. > > I cannot test the ld128 implementation, but I believe that it works > > (for some definition of work). > > > > I forgot to mention that the initial files under ld{80,128} where > created via > > svn cp src/k_exp.c ld80/k_cexpl.c > svn cp src/s_cexp.c ld80/s_cexpl.c > svn cp src/k_exp.c ld128/k_cexpl.c > svn cp src/s_cexp.c ld128/s_cexpl.c > > so the diff is not against an empty file; but against the > copied files. > Just realized that I forgot to update the cexp.3 manual page. -- Steve From owner-freebsd-numerics@freebsd.org Mon Feb 4 22:39:07 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 3C4BF14B519B for ; Mon, 4 Feb 2019 22:39:07 +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 BD08275DEF for ; Mon, 4 Feb 2019 22:39:05 +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 x14Md3Wl051260 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Mon, 4 Feb 2019 14:39:03 -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 x14Md3d5051259 for freebsd-numerics@freebsd.org; Mon, 4 Feb 2019 14:39:03 -0800 (PST) (envelope-from sgk) Date: Mon, 4 Feb 2019 14:39:03 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: Implement cexpl Message-ID: <20190204223903.GA51252@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20190203173056.GA42499@troutmask.apl.washington.edu> <20190203174259.GA42562@troutmask.apl.washington.edu> <20190204193247.GA50265@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20190204193247.GA50265@troutmask.apl.washington.edu> User-Agent: Mutt/1.11.2 (2019-01-07) X-Rspamd-Queue-Id: BD08275DEF X-Spamd-Bar: +++ Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [3.28 / 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.93)[0.932,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.67)[0.674,0]; 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]; DMARC_NA(0.00)[washington.edu]; NEURAL_SPAM_LONG(0.86)[0.859,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.12)[ip: (0.28), ipnet: 128.95.0.0/16(0.33), 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: Mon, 04 Feb 2019 22:39:07 -0000 On Mon, Feb 04, 2019 at 11:32:47AM -0800, Steve Kargl wrote: > > Just realized that I forgot to update the cexp.3 manual page. > Edits for cexp.3. Index: Makefile =================================================================== --- Makefile (revision 2143) +++ Makefile (working copy) @@ -198,7 +198,7 @@ MLINKS+=ccos.3 ccosf.3 ccos.3 csin.3 ccos.3 csinf.3 cc MLINKS+=ccosh.3 ccoshf.3 ccosh.3 csinh.3 ccosh.3 csinhf.3 \ ccosh.3 ctanh.3 ccosh.3 ctanhf.3 MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3 -MLINKS+=cexp.3 cexpf.3 +MLINKS+=cexp.3 cexpf.3 cexp.3 cexpl.3 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \ cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \ cimag.3 cproj.3 cimag.3 cprojf.3 cimag.3 cprojl.3 \ Index: man/cexp.3 =================================================================== --- man/cexp.3 (revision 2143) +++ man/cexp.3 (working copy) @@ -24,12 +24,13 @@ .\" .\" $FreeBSD: head/lib/msun/man/cexp.3 276293 2014-12-27 08:22:58Z joel $ .\" -.Dd March 6, 2011 +.Dd February 4, 2019 .Dt CEXP 3 .Os .Sh NAME .Nm cexp , -.Nm cexpf +.Nm cexpf , +.Nm cexpl .Nd complex exponential functions .Sh LIBRARY .Lb libm @@ -39,11 +40,14 @@ .Fn cexp "double complex z" .Ft float complex .Fn cexpf "float complex z" +.Ft long double complex +.Fn cexpl "long double complex z" .Sh DESCRIPTION The -.Fn cexp +.Fn cexp , +.Fn cexpf , and -.Fn cexpf +.Fn cexpl functions compute the complex exponential of .Fa z , also known as @@ -106,8 +110,9 @@ is not finite, the sign of the result is indeterminate .Xr math 3 .Sh STANDARDS The -.Fn cexp +.Fn cexp , +.Fn cexpf , and -.Fn cexpf +.Fn cexpl functions conform to .St -isoC-99 . -- Steve