From owner-freebsd-numerics@freebsd.org Wed Feb 27 02:01:46 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 94E81150BD45 for ; Wed, 27 Feb 2019 02:01:46 +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 421CC86A1E for ; Wed, 27 Feb 2019 02:01:45 +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 x1R21gEa075080 (version=TLSv1.3 cipher=TLS_AES_256_GCM_SHA384 bits=256 verify=NO) for ; Tue, 26 Feb 2019 18:01:42 -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 x1R21g2g075079 for freebsd-numerics@freebsd.org; Tue, 26 Feb 2019 18:01:42 -0800 (PST) (envelope-from sgk) Date: Tue, 26 Feb 2019 18:01:42 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: Implement cexpl Message-ID: <20190227020142.GA74833@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: 421CC86A1E X-Spamd-Bar: / Authentication-Results: mx1.freebsd.org X-Spamd-Result: default: False [0.18 / 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_HAM_LONG(-0.25)[-0.254,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)[]; 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_HAM_SHORT(-0.79)[-0.790,0]; NEURAL_SPAM_MEDIUM(0.48)[0.483,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.05)[ip: (0.11), ipnet: 128.95.0.0/16(0.16), asn: 73(0.06), 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: Wed, 27 Feb 2019 02:01:46 -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 have attached a new patch to https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=235413 This new patch incorporates changes to cexp and cexpf to compute the z = x + I y with x == 0 case first. This is in preparation for fixing constanting folding with GCC. The implementation for ld128/*_cexpl.c have not been compiled, and hence, are untested. I don't have ld128 hardware. The specific changes are * lib/msun/src/math_private.h: . Add an EXTRACT_LDBL80_WORDS2() macro to get access to the high and low word of a 64-bit significand as well as the expsign. . Add prototype for __ldexp_expl(). . Add prototype for __ldexp_cexpl(). * lib/msun/src/s_cexp.c: . A float.h to get LDBL_MANT_DIG. . Add c and s to declaration, and sort. . Move x = 0 case to be the first case tested. This is in preparation for fixing constanting folding in GCC. . Use sincos() instead of a call to sin() and to cos(). . A week_refrence for LDBL_MANT_DIG == 53. * lib/msun/src/s_cexpf.c: . Add c and s to declaration, and sort. . Move x = 0 case to be the first case tested. This is in preparation for fixing constanting folding in GCC. . Use sincosf() instead of a call to sinf() and to cosf(). * lib/msun/src/k_exp.c: . Add c and s to declaration, and sort. . Use sincos() instead of a call to sin() and to cos(). * lib/msun/src/k_expf.c . Add c and s to declaration, and sort. . Use sincosf() instead of a call to sinf() and to cosf(). * lib/msun/ld128/k_cexpl.c: . Copy src/k_exp.c to here. #if 0 ... #endif all code and have functions return NaN. These functions are currently unused, and await someone who cares. . Issue a compile-time warning about the sloppiness. * lib/msun/ld128/s_cexpl.c: . Copy src/s_cexp.c to here. . Convert "double complex" to "long double complex" without use of bit-twiddling. . Add compile-time warning about the sloppiness. . Add run-time warning about the sloppiness. * lib/msun/ld80/k_cexpl.c: . Copy src/k_exp.c to here. . Convert "double complex" to "long double complex". Use bit-twiddling. * lib/msun/ld80/s_cexpl.c: . Copy src/s_cexp.c to here. . Convert "double complex" to "long double complex". Use bit-twiddling where bits are grabbed with new EXTRACT_LDBL80_WORDS2() macro. * lib/msun/man/cexp.3: . Document the addtion of cexpl. * include/complex.h: . Add prototype for cexpl(). Index: lib/msun/src/math_private.h =================================================================== --- lib/msun/src/math_private.h (revision 344600) +++ lib/msun/src/math_private.h (working copy) @@ -243,6 +243,20 @@ } while (0) /* + * Get expsign and mantissa as 16 bit and 2 32-bit ints from an 80 bit long + * double. + */ + +#define EXTRACT_LDBL80_WORDS2(ix0,ix1,ix2,d) \ +do { \ + union IEEEl2bits ew_u; \ + ew_u.e = (d); \ + (ix0) = ew_u.xbits.expsign; \ + (ix1) = ew_u.bits.manh; \ + (ix2) = ew_u.bits.manl; \ +} while (0) + +/* * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit * long double. */ @@ -920,5 +934,9 @@ 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 #endif /* !_MATH_PRIVATE_H_ */ Index: lib/msun/src/s_cexp.c =================================================================== --- lib/msun/src/s_cexp.c (revision 344600) +++ lib/msun/src/s_cexp.c (working copy) @@ -30,6 +30,7 @@ __FBSDID("$FreeBSD$"); #include +#include #include #include "math_private.h" @@ -41,12 +42,19 @@ 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); y = cimag(z); + 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)); + } + EXTRACT_WORDS(hy, ly, y); hy &= 0x7fffffff; @@ -53,10 +61,6 @@ /* 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) - return (CMPLX(cos(y), sin(y))); if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { @@ -86,6 +90,11 @@ * - 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)); } } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(cexp, cexpl); +#endif Index: lib/msun/src/s_cexpf.c =================================================================== --- lib/msun/src/s_cexpf.c (revision 344600) +++ lib/msun/src/s_cexpf.c (working copy) @@ -41,12 +41,19 @@ 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); y = cimagf(z); + 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; @@ -53,10 +60,6 @@ /* 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) - return (CMPLXF(cosf(y), sinf(y))); if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { @@ -86,6 +89,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)); } } Index: lib/msun/src/k_exp.c =================================================================== --- lib/msun/src/k_exp.c (revision 344600) +++ lib/msun/src/k_exp.c (working copy) @@ -88,7 +88,7 @@ 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 @@ 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: lib/msun/src/k_expf.c =================================================================== --- lib/msun/src/k_expf.c (revision 344600) +++ lib/msun/src/k_expf.c (working copy) @@ -71,7 +71,7 @@ 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 @@ 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: lib/msun/ld128/k_cexpl.c =================================================================== --- lib/msun/ld128/k_cexpl.c (nonexistent) +++ lib/msun/ld128/k_cexpl.c (working copy) @@ -0,0 +1,126 @@ +/*- + * 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. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include + +#include "math.h" +#include "math_private.h" + +#warning These functions are 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 + * returned separately in 'expt'. + * + * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 + * Output: 2**1023 <= y < 2**1024 + */ +static long double +__frexp_expl(long double x, int *expt) +{ +#if 0 + double exp_x; + uint32_t hx; + + /* + * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to + * minimize |exp(kln2) - 2**k|. We also scale the exponent of + * 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); +#endif + return (x - x) / (x - x); +} + +/* + * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt. + * They are intended for large arguments (real part >= ln(DBL_MAX)) + * where care is needed to avoid overflow. + * + * The present implementation is narrowly tailored for our hyperbolic and + * exponential functions. We assume expt is small (0 or -1), and the caller + * has filtered out very large x, for which overflow would be inevitable. + */ + +long double +__ldexp_expl(long double x, int expt) +{ +#if 0 + double exp_x, scale; + int ex_expt; + + exp_x = __frexp_exp(x, &ex_expt); + expt += ex_expt; + INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); + return (exp_x * scale); +#endif + return (x - x) / (x - x); +} + +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; + + x = creal(z); + y = cimag(z); + exp_x = __frexp_exp(x, &ex_expt); + expt += ex_expt; + + /* + * Arrange so that scale1 * scale2 == 2**expt. We use this to + * compensate for scalbn being horrendously slow. + */ + half_expt = expt / 2; + INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); + half_expt = expt - half_expt; + INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); + + sincos(y, &s, &c); + return (CMPLX(c) * exp_x * scale1 * scale2, + s * exp_x * scale1 * scale2)); +#endif + return (x - x) / (x - x); +} Property changes on: lib/msun/ld128/k_cexpl.c ___________________________________________________________________ Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Index: lib/msun/ld128/s_cexpl.c =================================================================== --- lib/msun/ld128/s_cexpl.c (nonexistent) +++ lib/msun/ld128/s_cexpl.c (working copy) @@ -0,0 +1,76 @@ +/*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * + * Copyright (c) 2019 Steven G. Kargl + * 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. + */ + +#include +__FBSDID("$FreeBSD$"); + +#include +#include +#include + +#include "math_private.h" + +#warning cexpl is likely broken on ld128 architectures. +#warning Someone who cares needs to convert src/s_cexp.c. + +long double complex +cexpl(long double complex z) +{ + long double c, exp_x, s, x, y; + + x = creall(z); + y = cimagl(z); + + /* cexp(0 + I y) = cos(y) + I sin(y) */ + if (x == 0) { + sincosl(y, &s, &c); + return (CMPLXL(c, s)); + } + + /* 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 (CMPLXL(y - y, y - y)); + } else if (isinf(x) && copysignl(1.L, x) < 0) { + /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ + return (CMPLXL(0.0L, 0.0L)); + } else { + /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ + return (CMPLXL(x, y - y)); + } + } + + exp_x = expl(x); + sincosl(y, &s, &c); + return (CMPLXL(exp_x * c, exp_x * s)); +} +__warn_references("Precision of cexpl may have lower than expected precision"); Property changes on: lib/msun/ld128/s_cexpl.c ___________________________________________________________________ Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Index: lib/msun/ld80/k_cexpl.c =================================================================== --- lib/msun/ld80/k_cexpl.c (nonexistent) +++ lib/msun/ld80/k_cexpl.c (working copy) @@ -0,0 +1,118 @@ +/*- + * 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/k_exp.c converted to long double complex by Steven G. Kargl + */ + +#include +__FBSDID("$FreeBSD$"); + +#include + +#include "fpmath.h" +#include "math.h" +#include "math_private.h" + +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(LDBL_MAX) <= x < ln(2 * LDBL_MAX / LDBL_MIN_DENORM) ~= 22756.0 + * Output: 2**16383 <= y < 2**16384 + */ +static long double +__frexp_expl(long double x, int *expt) +{ + union IEEEl2bits exp_x; + + /* + * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to + * minimize |exp(kln2) - 2**k|. We also scale the exponent of + * exp_x to MAX_EXP so that the result can be multiplied by + * a tiny number without losing accuracy due to denormalization. + */ + 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); +} + +/* + * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt. + * They are intended for large arguments (real part >= ln(DBL_MAX)) + * where care is needed to avoid overflow. + * + * The present implementation is narrowly tailored for our hyperbolic and + * exponential functions. We assume expt is small (0 or -1), and the caller + * has filtered out very large x, for which overflow would be inevitable. + */ + +long double +__ldexp_expl(long double x, int expt) +{ + union IEEEl2bits scale; + long double exp_x; + int ex_expt; + + exp_x = __frexp_expl(x, &ex_expt); + expt += ex_expt; + scale.e = 1; + scale.bits.exp = 0x3fff + expt; + return (exp_x * scale.e); +} + +long double complex +__ldexp_cexpl(long double complex z, int expt) +{ + union IEEEl2bits scale1, scale2; + long double c, exp_x, s, x, y; + int ex_expt, half_expt; + + x = creall(z); + y = cimagl(z); + exp_x = __frexp_expl(x, &ex_expt); + expt += ex_expt; + + /* + * Arrange so that scale1 * scale2 == 2**expt. We use this to + * compensate for scalbn being horrendously slow. + */ + half_expt = expt / 2; + scale1.e = 1; + scale1.bits.exp = 0x3fff + half_expt; + half_expt = expt - half_expt; + scale2.e = 1; + scale2.bits.exp = 0x3fff + half_expt + 1; + + sincosl(y, &s, &c); + return (CMPLXL(c * exp_x * scale1.e * scale2.e, + s * exp_x * scale1.e * scale2.e)); +} Property changes on: lib/msun/ld80/k_cexpl.c ___________________________________________________________________ Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Index: lib/msun/ld80/s_cexpl.c =================================================================== --- lib/msun/ld80/s_cexpl.c (nonexistent) +++ lib/msun/ld80/s_cexpl.c (working copy) @@ -0,0 +1,105 @@ +/*- + * 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 + +#include "fpmath.h" +#include "math_private.h" + +static const uint32_t +exp_ovfl = 0xb17217f7, /* high bits of MAX_EXP * ln2 ~= 11356 */ +cexp_ovfl = 0xb1c6a857; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ + +long double complex +cexpl (long double complex z) +{ + long double c, exp_x, s, x, y; + uint32_t hwx, hwy, lwx, lwy; + uint16_t hx, hy; + + ENTERI(z); + + x = creall(z); + y = cimagl(z); + + EXTRACT_LDBL80_WORDS2(hx, hwx, lwx, x); + EXTRACT_LDBL80_WORDS2(hy, hwy, lwy, y); + + /* cexp(0 + I y) = cos(y) + I sin(y) */ + if ((hwx | lwx | (hx & 0x7fff)) == 0) { + sincosl(y, &s, &c); + RETURNI(CMPLXL(c, s)); + } + + /* cexp(x + I 0) = exp(x) + I 0 */ + if ((hwy | lwy | (hy & 0x7fff)) == 0) + RETURNI(CMPLXL(expl(x), y)); + + if (hy >= 0x7fff) { + if (hx < 0x7fff || + (hx == 0x7fff && (hwx & 0x7fffffff) != 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.0L, 0.0L)); + } else { + /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ + RETURNI(CMPLXL(x, y - y)); + } + } + + if (hx >= 0x400c && hx <= 0x400d) { + /* + * x is between 11356 and 22755, so we must scale to avoid + * overflow in exp(x). + */ + RETURNI(__ldexp_cexpl(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 = expl(x); + sincosl(y, &s, &c); + RETURNI(CMPLXL(exp_x * c, exp_x * s)); + } +} Property changes on: lib/msun/ld80/s_cexpl.c ___________________________________________________________________ Added: svn:eol-style ## -0,0 +1 ## +native \ No newline at end of property Added: svn:keywords ## -0,0 +1 ## +FreeBSD=%H \ No newline at end of property Added: svn:mime-type ## -0,0 +1 ## +text/plain \ No newline at end of property Index: lib/msun/man/cexp.3 =================================================================== --- lib/msun/man/cexp.3 (revision 344600) +++ lib/msun/man/cexp.3 (working copy) @@ -24,12 +24,13 @@ .\" .\" $FreeBSD$ .\" -.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 @@ .Xr math 3 .Sh STANDARDS The -.Fn cexp +.Fn cexp , +.Fn cexpf , and -.Fn cexpf +.Fn cexpl functions conform to .St -isoC-99 . Index: include/complex.h =================================================================== --- include/complex.h (revision 344600) +++ include/complex.h (working copy) @@ -98,6 +98,8 @@ 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;