Attachment 202398 Details for Bug 235413 – New patch with cexpl implementation included

[patch] New patch with cexpl implementation included

cexpl.diff (text/plain), 21.83 KB, created by Steve Kargl on 2019-02-27 01:56:55 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Steve Kargl

Created: 2019-02-27 01:56:55 UTC

Size: 21.83 KB

patch

obsolete

>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@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@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@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@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;

Actions: View | Diff

Attachments on bug 235413: 201620 | 202398 | 208609