Implement cexpl
Steve Kargl
sgk at troutmask.apl.washington.edu
Wed Feb 27 02:01:46 UTC 2019
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 <complex.h>
+#include <float.h>
#include <math.h>
#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 <das at FreeBSD.ORG>
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+
+#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 <kargl at FreeBSD.ORG>
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <float.h>
+#include <math.h>
+
+#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 <das at FreeBSD.ORG>
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+
+#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 <das at FreeBSD.ORG>
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+#include <math.h>
+
+#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;
More information about the freebsd-numerics
mailing list