Attachment 208609 Details for Bug 235413 – Updated patch

[patch] Updated patch

cexpl_new.diff (text/plain), 24.30 KB, created by Steve Kargl on 2019-10-26 16:32:28 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Steve Kargl

Created: 2019-10-26 16:32:28 UTC

Size: 24.30 KB

patch

obsolete

>Index: include/complex.h
>===================================================================
>--- include/complex.h	(revision 354057)
>+++ 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;
>Index: lib/msun/Makefile
>===================================================================
>--- lib/msun/Makefile	(revision 354057)
>+++ lib/msun/Makefile	(working copy)
>@@ -101,7 +101,7 @@
> 	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 \
>-	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c \
>+	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cexpl.c \
> 	s_clogl.c s_cosl.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 \
>@@ -173,7 +173,7 @@
> 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: lib/msun/Symbol.map
>===================================================================
>--- lib/msun/Symbol.map	(revision 354057)
>+++ lib/msun/Symbol.map	(working copy)
>@@ -304,3 +304,8 @@
> 	sincosf;
> 	sincosl;
> };
>+
>+/* First added in 13.0-CURRENT */
>+FBSD_1.6 {
>+	cexpl;
>+};
>Index: lib/msun/src/math_private.h
>===================================================================
>--- lib/msun/src/math_private.h	(revision 354057)
>+++ lib/msun/src/math_private.h	(working copy)
>@@ -335,18 +334,19 @@
> 
> /* Support switching the mode to FP_PE if necessary. */
> #if defined(__i386__) && !defined(NO_FPSETPREC)
>-#define	ENTERI() ENTERIT(long double)
>-#define	ENTERIT(returntype)			\
>-	returntype __retval;			\
>+#define	ENTERI()				\
> 	fp_prec_t __oprec;			\
> 						\
> 	if ((__oprec = fpgetprec()) != FP_PE)	\
> 		fpsetprec(FP_PE)
>-#define	RETURNI(x) do {				\
>-	__retval = (x);				\
>+#define	LEAVEI()				\
> 	if (__oprec != FP_PE)			\
>-		fpsetprec(__oprec);		\
>-	RETURNF(__retval);			\
>+		fpsetprec(__oprec)
>+#define	RETURNI(expr) do {                      \
>+	__typeof(expr) __retval = (expr);	\
>+						\
>+	LEAVEI();                               \
>+	RETURNF(__retval);                      \
> } while (0)
> #define	ENTERV()				\
> 	fp_prec_t __oprec;			\
>@@ -356,12 +356,12 @@
> #define	RETURNV() do {				\
> 	if (__oprec != FP_PE)			\
> 		fpsetprec(__oprec);		\
>-	return;			\
>+	return;					\
> } while (0)
> #else
> #define	ENTERI()
>-#define	ENTERIT(x)
>-#define	RETURNI(x)	RETURNF(x)
>+#define	LEAVEI()
>+#define	RETURNI(expr)	RETURNF(expr)
> #define	ENTERV()
> #define	RETURNV()	return
> #endif
>@@ -920,5 +920,6 @@
> 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);
> 
> #endif /* !_MATH_PRIVATE_H_ */
>Index: lib/msun/man/cexp.3
>===================================================================
>--- lib/msun/man/cexp.3	(revision 354057)
>+++ 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: lib/msun/src/k_exp.c
>===================================================================
>--- lib/msun/src/k_exp.c	(revision 354057)
>+++ 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 354057)
>+++ 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/src/s_cexp.c
>===================================================================
>--- lib/msun/src/s_cexp.c	(revision 354057)
>+++ 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,7 +42,7 @@
> double complex
> cexp(double complex z)
> {
>-	double x, y, exp_x;
>+	double c, exp_x, s, x, y;
> 	uint32_t hx, hy, lx, ly;
> 
> 	x = creal(z);
>@@ -55,8 +56,10 @@
> 		return (CMPLX(exp(x), y));
> 	EXTRACT_WORDS(hx, lx, x);
> 	/* cexp(0 + I y) = cos(y) + I sin(y) */
>-	if (((hx & 0x7fffffff) | lx) == 0)
>-		return (CMPLX(cos(y), sin(y)));
>+	if (((hx & 0x7fffffff) | lx) == 0) {
>+		sincos(y, &s, &c);
>+		return (CMPLX(c, s));
>+	}
> 
> 	if (hy >= 0x7ff00000) {
> 		if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
>@@ -86,6 +89,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 354057)
>+++ lib/msun/src/s_cexpf.c	(working copy)
>@@ -41,7 +41,7 @@
> float complex
> cexpf(float complex z)
> {
>-	float x, y, exp_x;
>+	float c, exp_x, s, x, y;
> 	uint32_t hx, hy;
> 
> 	x = crealf(z);
>@@ -55,8 +55,10 @@
> 		return (CMPLXF(expf(x), y));
> 	GET_FLOAT_WORD(hx, x);
> 	/* cexp(0 + I y) = cos(y) + I sin(y) */
>-	if ((hx & 0x7fffffff) == 0)
>-		return (CMPLXF(cosf(y), sinf(y)));
>+	if ((hx & 0x7fffffff) == 0) {
>+		sincosf(y, &s, &c);
>+		return (CMPLXF(c, s));
>+	}
> 
> 	if (hy >= 0x7f800000) {
> 		if ((hx & 0x7fffffff) != 0x7f800000) {
>@@ -86,6 +88,7 @@
> 		 *  -  x = NaN (spurious inexact exception from y)
> 		 */
> 		exp_x = expf(x);
>-		return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
>+		sincosf(y, &s, &c);
>+		return (CMPLXF(exp_x * c, exp_x * s));
> 	}
> }
>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));
>+}
>
>Index: lib/msun/ld80/k_expl.h
>===================================================================
>--- lib/msun/ld80/k_expl.h	(revision 354057)
>+++ lib/msun/ld80/k_expl.h	(working copy)
>@@ -277,7 +277,7 @@
> static inline long double complex
> __ldexp_cexpl(long double complex z, int expt)
> {
>-	long double exp_x, hi, lo;
>+	long double c, exp_x, hi, lo, s;
> 	long double x, y, scale1, scale2;
> 	int half_expt, k;
> 
>@@ -285,7 +285,7 @@
> 	y = cimagl(z);
> 	__k_expl(x, &hi, &lo, &k);
> 
>-	exp_x = (lo + hi) * 0x1p16382;
>+	exp_x = (lo + hi) * 0x1p16382L;
> 	expt += k - 16382;
> 
> 	scale1 = 1;
>@@ -292,9 +292,10 @@
> 	half_expt = expt / 2;
> 	SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
> 	scale2 = 1;
>-	SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
>+	SET_LDBL_EXPSIGN(scale2, BIAS + expt - half_expt);
> 
>-	return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
>-	    sinl(y) * exp_x * scale1 * scale2));
>+	sincosl(y, &s, &c);
>+	return (CMPLXL(c * exp_x * scale1 * scale2,
>+	    s * exp_x * scale1 * scale2));
> }
> #endif /* _COMPLEX_H */
>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,107 @@
>+/*-
>+ * 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 "fpmath.h"
>+#include "math.h"
>+#include "math_private.h"
>+#include "k_expl.h"
>+
>+long double complex
>+cexpl (long double complex z)
>+{
>+	long double c, exp_x, s, x, y;
>+	uint64_t lx, ly;
>+	uint16_t hx, hy;
>+
>+	ENTERI();
>+
>+	x = creall(z);
>+	y = cimagl(z);
>+
>+	EXTRACT_LDBL80_WORDS(hy, ly, y);
>+	hy &= 0x7fff;
>+
>+	/* cexp(x + I 0) = exp(x) + I 0 */
>+	if ((hy | ly) == 0)
>+		RETURNI(CMPLXL(expl(x), y));
>+	EXTRACT_LDBL80_WORDS(hx, lx, x);
>+	/* cexp(0 + I y) = cos(y) + I sin(y) */
>+	if (((hx & 0x7fff) | lx) == 0) {
>+		sincosl(y, &s, &c);
>+		RETURNI(CMPLXL(c, s));
>+	}
>+
>+	if (hy >= 0x7fff) {
>+		if ((hx & 0x7fff) < 0x7fff || ((hx & 0x7fff) == 0x7fff &&
>+		    (lx & 0x7fffffffffffffffULL) != 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.0, 0.0));
>+		} else {
>+			/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
>+			RETURNI(CMPLXL(x, y - y));
>+		}
>+	}
>+
>+	/*
>+	 *  exp_ovfl = 11356.5234062941439497
>+	 * cexp_ovfl = 22755.3287906024445633
>+	 */
>+	if ((hx == 0x400c && lx > 0xb17217f7d1cf79acULL) ||
>+	    (hx == 0x400d && lx < 0xb1c6a8573de9768cULL)) {
>+		/*
>+		 * x is between exp_ovfl and cexp_ovfl, 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));
>+	}
>+}
>
>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);
>+}
>
>Index: lib/msun/ld128/k_expl.h
>===================================================================
>--- lib/msun/ld128/k_expl.h	(revision 354057)
>+++ lib/msun/ld128/k_expl.h	(working copy)
>@@ -299,7 +299,7 @@
> static inline long double complex
> __ldexp_cexpl(long double complex z, int expt)
> {
>-	long double exp_x, hi, lo;
>+	long double c, exp_x, hi, lo, s;
> 	long double x, y, scale1, scale2;
> 	int half_expt, k;
> 
>@@ -307,7 +307,7 @@
> 	y = cimagl(z);
> 	__k_expl(x, &hi, &lo, &k);
> 
>-	exp_x = (lo + hi) * 0x1p16382;
>+	exp_x = (lo + hi) * 0x1p16382L;
> 	expt += k - 16382;
> 
> 	scale1 = 1;
>@@ -314,9 +314,10 @@
> 	half_expt = expt / 2;
> 	SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
> 	scale2 = 1;
>-	SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
>+	SET_LDBL_EXPSIGN(scale2, BIAS + expt - half_expt);
> 
>-	return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
>-	    sinl(y) * exp_x * scale1 * scale2));
>+	sincosl(y, &s, &c);
>+	return (CMPLXL(c * exp_x * scale1 * scale2,
>+	    s * exp_x * scale1 * scale2));
> }
> #endif /* _COMPLEX_H */
>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,92 @@
>+/*-
>+ * 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() should be converted to use bits likeo src/s_cexp.c.
>+
>+static const long double
>+cexp_ovfl = 2.27892930024498818830197576893019292e+04L,
>+exp_ovfl = 1.13565234062941439494919310779707649e+04L;
>+
>+long double complex
>+cexpl(long double complex z)
>+{
>+	long double c, exp_x, s, x, y;
>+
>+	x = creall(z);
>+	y = cimagl(z);
>+
>+	/* cexp(x + I 0) = exp(x) + I 0 */
>+	if (y == 0)
>+		return (CMPLXL(expl(x), y));
>+	/* cexp(0 + I y) = cos(y) + I sin(y) */
>+	if (x == 0) {
>+		sincosl(y, &s, &c);
>+		return (CMPLXL(c, s));
>+	}
>+
>+	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.0, 0.0));
>+		} else {
>+			/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
>+			return (CMPLXL(x, y - y));
>+		}
>+	}
>+
>+	if (x > exp_ovfl && x < cexp_ovfl) {
>+		/*
>+		 * x is between exp_ovfl and cexp_ovfl, 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);
>+		return (CMPLXL(exp_x * c, exp_x * s));
>+	}
>+}

Actions: View | Diff

Attachments on bug 235413: 201620 | 202398 | 208609