Date: Sun, 3 Feb 2019 09:30:56 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: Implement cexpl Message-ID: <20190203173056.GA42499@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
The following patch implements cexpl() for ld80 and ld128 architectures.
The kernel files for large x are named ld{80,128}/k_cexpl.c to avoid
confusion with the ld{80,128}/k_expl.h files.
I do not have access to ld128, so the ld128 does not use bit twiddling.
I cannot test the ld128 implementation, but I believe that it works
(for some definition of work).
The large patch does not contain the needed change to include/complex.h.
Index: include/complex.h
===================================================================
--- include/complex.h (revision 343477)
+++ include/complex.h (working copy)
@@ -98,6 +98,8 @@ double complex ccosh(double complex);
float complex ccoshf(float complex);
double complex cexp(double complex);
float complex cexpf(float complex);
+long double complex
+ cexpl(long double complex);
double cimag(double complex) __pure2;
float cimagf(float complex) __pure2;
long double cimagl(long double complex) __pure2;
Here's a svn logfile.
* Makefile:
. Add k_cexpl.c and s_cexpl.c to build.
* src/math_private.h:
. Add ENTERC and RETURNC to allow toggling of FP_PE in i686-class
hardware for long double complex.
. Add prototype for __ldexp_cexpl().
* src/s_cexp.c:
. Include float.h to get LDBL_MANT_DIG so that a weak reference can be
add for LDBL_MANT_DIG == 53 targets.
. Flip the y == 0 and x == 0 cases. This is in anticipation of allowing
GCC to transform CMPLX(cos(x), sin(x)) into cexp(CMPLX(0,x)), which
can then be mapped to the sincos() builtin.
* src/s_cexpf.c:
. Flip the y == 0 and x == 0 cases. This is in anticipation of allowing
GCC to transform CMPLX(cos(x), sin(x)) into cexp(CMPLX(0,x)), which
can then be mapped to the sincos() builtin.
* src/k_exp.c:
. Declare c and s, and re-order declarations.
. Use sincos().
* src/k_expf.c:
. Declare c and s, and re-order declarations.
. Use sincosf().
* ld80/k_cexpl.c:
. Implementation of long double complex kernels for cexpl() for x large.
. Conversion of src/k_cexp.c from double complex to long double complex.
* ld80/s_cexpl.c:
. Implementation of long double complex cexpl().
. Conversion of src/k_cexp.c from double complex to long double complex.
* ld128/k_cexpl.c:
. FIXME: I do not have access to ld128 hardware. Someone who cares about
ld128 needs to fix this implementation of long double complex kernels
for cexpl() for x large.
. The functions defined here are currently unused.
* ld128/s_cexpl.c:
. FIXME: I do not have access to ld128 hardware. Someone who cares about
ld128 needs to fix this implementation of long double complex cexpl()
to use bit twiddling. It currently uses isinf() and isnan() where
needed and floating point comparisons.
And now the diff.
Index: Makefile
===================================================================
--- Makefile (revision 2141)
+++ Makefile (working copy)
@@ -133,7 +133,9 @@ COMMON_SRCS+= catrigl.c \
e_lgammal.c e_lgammal_r.c e_powl.c \
e_remainderl.c e_sinhl.c e_sqrtl.c \
invtrig.c k_cosl.c k_sinl.c k_tanl.c \
+ k_cexpl.c \
s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c \
+ s_cexpl.c \
s_clogl.c s_cosl.c s_cospil.c s_cprojl.c \
s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
Index: src/math_private.h
===================================================================
--- src/math_private.h (revision 2141)
+++ src/math_private.h (working copy)
@@ -357,9 +357,24 @@ do { \
fpsetprec(__oprec); \
return; \
} while (0)
+#define ENTERC() ENTERCT(long double complex)
+#define ENTERCT(returntype) \
+ returntype __retval; \
+ fp_prec_t __oprec; \
+ \
+ if ((__oprec = fpgetprec()) != FP_PE) \
+ fpsetprec(FP_PE)
+#define RETURNC(x) do { \
+ __retval = (x); \
+ if (__oprec != FP_PE) \
+ fpsetprec(__oprec); \
+ RETURNF(__retval); \
+} while (0)
#else
+#define ENTERC(x)
#define ENTERI(x)
#define ENTERIT(x)
+#define RETURNC(x) RETURNF(x)
#define RETURNI(x) RETURNF(x)
#define ENTERV()
#define RETURNV() return
@@ -919,6 +934,10 @@ float complex __ldexp_cexpf(float complex,int);
long double __kernel_sinl(long double, long double, int);
long double __kernel_cosl(long double, long double);
long double __kernel_tanl(long double, long double, int);
+long double __ldexp_expl(long double, int);
+#ifdef _COMPLEX_H
+long double complex __ldexp_cexpl(long double complex, int);
+#endif
#include "math_sgk.h"
Index: src/s_cexp.c
===================================================================
--- src/s_cexp.c (revision 2141)
+++ src/s_cexp.c (working copy)
@@ -30,6 +30,7 @@
__FBSDID("$FreeBSD: head/lib/msun/src/s_cexp.c 326219 2017-11-26 02:00:33Z pfg $");
#include <complex.h>
+#include <float.h>
#include <math.h>
#include "math_private.h"
@@ -47,12 +48,6 @@ cexp(double complex z)
x = creal(z);
y = cimag(z);
- EXTRACT_WORDS(hy, ly, y);
- hy &= 0x7fffffff;
-
- /* cexp(x + I 0) = exp(x) + I 0 */
- if ((hy | ly) == 0)
- return (CMPLX(exp(x), y));
EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if (((hx & 0x7fffffff) | lx) == 0) {
@@ -60,6 +55,13 @@ cexp(double complex z)
return (CMPLX(c, s));
}
+ EXTRACT_WORDS(hy, ly, y);
+ hy &= 0x7fffffff;
+
+ /* cexp(x + I 0) = exp(x) + I 0 */
+ if ((hy | ly) == 0)
+ return (CMPLX(exp(x), y));
+
if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
@@ -92,3 +94,7 @@ cexp(double complex z)
return (CMPLX(exp_x * c, exp_x * s));
}
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cexp, cexpl);
+#endif
Index: src/s_cexpf.c
===================================================================
--- src/s_cexpf.c (revision 2141)
+++ src/s_cexpf.c (working copy)
@@ -47,18 +47,19 @@ cexpf(float complex z)
x = crealf(z);
y = cimagf(z);
- GET_FLOAT_WORD(hy, y);
- hy &= 0x7fffffff;
-
- /* cexp(x + I 0) = exp(x) + I 0 */
- if (hy == 0)
- return (CMPLXF(expf(x), y));
GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if ((hx & 0x7fffffff) == 0) {
sincosf(y, &s, &c);
return (CMPLXF(c, s));
}
+
+ GET_FLOAT_WORD(hy, y);
+ hy &= 0x7fffffff;
+
+ /* cexp(x + I 0) = exp(x) + I 0 */
+ if (hy == 0)
+ return (CMPLXF(expf(x), y));
if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) {
Index: src/k_exp.c
===================================================================
--- src/k_exp.c (revision 2141)
+++ src/k_exp.c (working copy)
@@ -88,7 +88,7 @@ __ldexp_exp(double x, int expt)
double complex
__ldexp_cexp(double complex z, int expt)
{
- double x, y, exp_x, scale1, scale2;
+ double c, exp_x, s, scale1, scale2, x, y;
int ex_expt, half_expt;
x = creal(z);
@@ -105,6 +105,7 @@ __ldexp_cexp(double complex z, int expt)
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
- return (CMPLX(cos(y) * exp_x * scale1 * scale2,
- sin(y) * exp_x * scale1 * scale2));
+ sincos(y, &s, &c);
+ return (CMPLX(c * exp_x * scale1 * scale2,
+ s * exp_x * scale1 * scale2));
}
Index: src/k_expf.c
===================================================================
--- src/k_expf.c (revision 2141)
+++ src/k_expf.c (working copy)
@@ -71,7 +71,7 @@ __ldexp_expf(float x, int expt)
float complex
__ldexp_cexpf(float complex z, int expt)
{
- float x, y, exp_x, scale1, scale2;
+ float c, exp_x, s, scale1, scale2, x, y;
int ex_expt, half_expt;
x = crealf(z);
@@ -84,6 +84,7 @@ __ldexp_cexpf(float complex z, int expt)
half_expt = expt - half_expt;
SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
- return (CMPLXF(cosf(y) * exp_x * scale1 * scale2,
- sinf(y) * exp_x * scale1 * scale2));
+ sincosf(y, &s, &c);
+ return (CMPLXF(c * exp_x * scale1 * scale2,
+ s * exp_x * scale1 * scale2));
}
Index: ld80/k_cexpl.c
===================================================================
--- ld80/k_cexpl.c (revision 2141)
+++ ld80/k_cexpl.c (working copy)
@@ -24,31 +24,34 @@
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * src/k_exp.c converted to long double complex by Steven G. Kargl
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/k_exp.c 326219 2017-11-26 02:00:33Z pfg $");
+__FBSDID("$FreeBSD$");
#include <complex.h>
+#include "fpmath.h"
#include "math.h"
#include "math_private.h"
-static const uint32_t k = 1799; /* constant for reduction */
-static const double kln2 = 1246.97177782734161156; /* k * ln2 */
-
+static const uint32_t k = 22956; /* constant for reduction */
+static const union IEEEl2bits
+kln2u = LD80C(0xf89f8bf509c862ca, 13, 1.59118866769341045231e+04L);
+#define kln2 (kln2u.e)
/*
* Compute exp(x), scaled to avoid spurious overflow. An exponent is
* returned separately in 'expt'.
*
- * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
- * Output: 2**1023 <= y < 2**1024
+ * Input: ln(LDBL_MAX) <= x < ln(2 * LDBL_MAX / LDBL_MIN_DENORM) ~= 22756.0
+ * Output: 2**16383 <= y < 2**16384
*/
-static double
-__frexp_exp(double x, int *expt)
+static long double
+__frexp_expl(long double x, int *expt)
{
- double exp_x;
- uint32_t hx;
+ union IEEEl2bits exp_x;
/*
* We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
@@ -56,11 +59,10 @@ __frexp_exp(double x, int *expt)
* exp_x to MAX_EXP so that the result can be multiplied by
* a tiny number without losing accuracy due to denormalization.
*/
- exp_x = exp(x - kln2);
- GET_HIGH_WORD(hx, exp_x);
- *expt = (hx >> 20) - (0x3ff + 1023) + k;
- SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
- return (exp_x);
+ exp_x.e = expl(x - kln2);
+ *expt = exp_x.bits.exp - (0x3fff + 16383) + k - 1;
+ exp_x.bits.exp = 0x3fff + 16383;
+ return (exp_x.e);
}
/*
@@ -73,27 +75,30 @@ __frexp_exp(double x, int *expt)
* has filtered out very large x, for which overflow would be inevitable.
*/
-double
-__ldexp_exp(double x, int expt)
+long double
+__ldexp_expl(long double x, int expt)
{
- double exp_x, scale;
+ union IEEEl2bits scale;
+ long double exp_x;
int ex_expt;
- exp_x = __frexp_exp(x, &ex_expt);
+ exp_x = __frexp_expl(x, &ex_expt);
expt += ex_expt;
- INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);
- return (exp_x * scale);
+ scale.e = 1;
+ scale.bits.exp = 0x3fff + expt;
+ return (exp_x * scale.e);
}
-double complex
-__ldexp_cexp(double complex z, int expt)
+long double complex
+__ldexp_cexpl(long double complex z, int expt)
{
- double x, y, exp_x, scale1, scale2;
+ union IEEEl2bits scale1, scale2;
+ long double c, exp_x, s, x, y;
int ex_expt, half_expt;
- x = creal(z);
- y = cimag(z);
- exp_x = __frexp_exp(x, &ex_expt);
+ x = creall(z);
+ y = cimagl(z);
+ exp_x = __frexp_expl(x, &ex_expt);
expt += ex_expt;
/*
@@ -101,10 +106,13 @@ __ldexp_cexp(double complex z, int expt)
* compensate for scalbn being horrendously slow.
*/
half_expt = expt / 2;
- INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0);
+ scale1.e = 1;
+ scale1.bits.exp = 0x3fff + half_expt;
half_expt = expt - half_expt;
- INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
+ scale2.e = 1;
+ scale2.bits.exp = 0x3fff + half_expt + 1;
- return (CMPLX(cos(y) * exp_x * scale1 * scale2,
- sin(y) * exp_x * scale1 * scale2));
+ sincosl(y, &s, &c);
+ return (CMPLXL(c * exp_x * scale1.e * scale2.e,
+ s * exp_x * scale1.e * scale2.e));
}
Index: ld80/s_cexpl.c
===================================================================
--- ld80/s_cexpl.c (revision 2141)
+++ ld80/s_cexpl.c (working copy)
@@ -24,61 +24,72 @@
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
+ *
+ * src/s_cexp.c converted to long double complex by Steven G. Kargl
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_cexp.c 326219 2017-11-26 02:00:33Z pfg $");
+__FBSDID("$FreeBSD$");
#include <complex.h>
+#include <float.h>
#include <math.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
#include "math_private.h"
static const uint32_t
-exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */
-cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
+exp_ovfl = 0xb17217f7, /* high bits of MAX_EXP * ln2 ~= 11356 */
+cexp_ovfl = 0xb1c6a857; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
-double complex
-cexp(double complex z)
+long double complex
+cexpl (long double complex z)
{
- double c, exp_x, s, x, y;
- uint32_t hx, hy, lx, ly;
+ union IEEEl2bits ux, uy;
+ long double c, exp_x, s, x, y;
+ uint16_t hx, hy;
- x = creal(z);
- y = cimag(z);
+ ENTERC();
- EXTRACT_WORDS(hy, ly, y);
- hy &= 0x7fffffff;
+ x = ux.e = creall(z);
+ y = uy.e = cimagl(z);
- /* cexp(x + I 0) = exp(x) + I 0 */
- if ((hy | ly) == 0)
- return (CMPLX(exp(x), y));
- EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
- if (((hx & 0x7fffffff) | lx) == 0) {
- sincos(y, &s, &c);
- return (CMPLX(c, s));
+ hx = ux.xbits.expsign;
+ if ((ux.bits.manh | ux.bits.manl | (hx & 0x7fff)) == 0) {
+ sincosl(y, &s, &c);
+ RETURNC(CMPLXL(c, s));
}
- if (hy >= 0x7ff00000) {
- if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
+ /* cexp(x + I 0) = exp(x) + I 0 */
+ hy = uy.xbits.expsign;
+ if ((uy.bits.manh | uy.bits.manl | (hy & 0x7fff)) == 0)
+ RETURNC(CMPLXL(expl(x), y));
+
+ if (hy >= 0x7fff) {
+ if (hx < 0x7fff ||
+ (hx == 0x7fff && (ux.bits.manh & 0x7fffffff) != 0)) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
- return (CMPLX(y - y, y - y));
- } else if (hx & 0x80000000) {
+ RETURNC(CMPLXL(y - y, y - y));
+ } else if (hx & 0x8000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
- return (CMPLX(0.0, 0.0));
+ RETURNC(CMPLXL(0.0L, 0.0L));
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
- return (CMPLX(x, y - y));
+ RETURNC(CMPLXL(x, y - y));
}
}
- if (hx >= exp_ovfl && hx <= cexp_ovfl) {
+ if (hx >= 0x400c && hx <= 0x400d) {
/*
- * x is between 709.7 and 1454.3, so we must scale to avoid
+ * x is between 11356 and 22755, so we must scale to avoid
* overflow in exp(x).
*/
- return (__ldexp_cexp(z, 0));
+ RETURNC(__ldexp_cexpl(z, 0));
} else {
/*
* Cases covered here:
@@ -87,8 +98,8 @@ cexp(double complex z)
* - x = +-Inf (generated by exp())
* - x = NaN (spurious inexact exception from y)
*/
- exp_x = exp(x);
- sincos(y, &s, &c);
- return (CMPLX(exp_x * c, exp_x * s));
+ exp_x = expl(x);
+ sincosl(y, &s, &c);
+ RETURNC(CMPLXL(exp_x * c, exp_x * s));
}
}
Index: ld128/k_cexpl.c
===================================================================
--- ld128/k_cexpl.c (revision 2141)
+++ ld128/k_cexpl.c (working copy)
@@ -34,8 +34,14 @@ __FBSDID("$FreeBSD: head/lib/msun/src/k_exp.c 326219 2
#include "math.h"
#include "math_private.h"
+#warning libm is broken on ld128 architectures
+#warning functions defined here are currently unused
+#warning someone who cares needs to convert src/k_exp.c
+
+#if 0
static const uint32_t k = 1799; /* constant for reduction */
static const double kln2 = 1246.97177782734161156; /* k * ln2 */
+#endif
/*
* Compute exp(x), scaled to avoid spurious overflow. An exponent is
@@ -44,9 +50,10 @@ static const double kln2 = 1246.97177782734161156; /*
* Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
* Output: 2**1023 <= y < 2**1024
*/
-static double
-__frexp_exp(double x, int *expt)
+static long double
+__frexp_expl(long double x, int *expt)
{
+#if 0
double exp_x;
uint32_t hx;
@@ -61,6 +68,8 @@ __frexp_exp(double x, int *expt)
*expt = (hx >> 20) - (0x3ff + 1023) + k;
SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
return (exp_x);
+#endif
+ return x * expt;
}
/*
@@ -73,9 +82,10 @@ __frexp_exp(double x, int *expt)
* has filtered out very large x, for which overflow would be inevitable.
*/
-double
-__ldexp_exp(double x, int expt)
+long double
+__ldexp_expl(long double x, int expt)
{
+#if 0
double exp_x, scale;
int ex_expt;
@@ -83,11 +93,14 @@ __ldexp_exp(double x, int expt)
expt += ex_expt;
INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);
return (exp_x * scale);
+#endif
+ return x * expt;
}
-double complex
-__ldexp_cexp(double complex z, int expt)
+long double complex
+__ldexp_cexpl(long double complex z, int expt)
{
+#if 0
double x, y, exp_x, scale1, scale2;
int ex_expt, half_expt;
@@ -105,6 +118,9 @@ __ldexp_cexp(double complex z, int expt)
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
- return (CMPLX(cos(y) * exp_x * scale1 * scale2,
- sin(y) * exp_x * scale1 * scale2));
+ sincos(y, &s, &c);
+ return (CMPLX(c) * exp_x * scale1 * scale2,
+ s * exp_x * scale1 * scale2));
+#endif
+ return z * expt;
}
Index: ld128/s_cexpl.c
===================================================================
--- ld128/s_cexpl.c (revision 2141)
+++ ld128/s_cexpl.c (working copy)
@@ -1,7 +1,7 @@
/*-
* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
*
- * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * Copyright (c) 2019 Steven G. Kargl <kargl@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
@@ -27,68 +27,49 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_cexp.c 326219 2017-11-26 02:00:33Z pfg $");
+__FBSDID("$FreeBSD$");
#include <complex.h>
+#include <float.h>
#include <math.h>
#include "math_private.h"
-static const uint32_t
-exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */
-cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
+#warning libm is broken on ld128 architectures
+#warning someone who cares needs to convert src/s_cexp.c
-double complex
-cexp(double complex z)
+long double complex
+cexpl(long double complex z)
{
- double c, exp_x, s, x, y;
- uint32_t hx, hy, lx, ly;
+ long double c, exp_x, s, x, y;
- x = creal(z);
- y = cimag(z);
+ x = creall(z);
+ y = cimagl(z);
- EXTRACT_WORDS(hy, ly, y);
- hy &= 0x7fffffff;
-
- /* cexp(x + I 0) = exp(x) + I 0 */
- if ((hy | ly) == 0)
- return (CMPLX(exp(x), y));
- EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
- if (((hx & 0x7fffffff) | lx) == 0) {
- sincos(y, &s, &c);
- return (CMPLX(c, s));
+ if (x == 0) {
+ sincosl(y, &s, &c);
+ return (CMPLXL(c, s));
}
- if (hy >= 0x7ff00000) {
- if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
+ /* cexp(x + I 0) = exp(x) + I 0 */
+ if (y == 0)
+ return (CMPLXL(expl(x), y));
+
+ if (!isfinite(y)) {
+ if (isfinite(x) || isnan(x)) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
- return (CMPLX(y - y, y - y));
- } else if (hx & 0x80000000) {
+ return (CMPLXL(y - y, y - y));
+ } else if (isinf(x) && copysignl(1.L, x) < 0) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
- return (CMPLX(0.0, 0.0));
+ return (CMPLXL(0.0, 0.0));
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
- return (CMPLX(x, y - y));
+ return (CMPLXL(x, y - y));
}
}
- if (hx >= exp_ovfl && hx <= cexp_ovfl) {
- /*
- * x is between 709.7 and 1454.3, so we must scale to avoid
- * overflow in exp(x).
- */
- return (__ldexp_cexp(z, 0));
- } else {
- /*
- * Cases covered here:
- * - x < exp_ovfl and exp(x) won't overflow (common case)
- * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
- * - x = +-Inf (generated by exp())
- * - x = NaN (spurious inexact exception from y)
- */
- exp_x = exp(x);
- sincos(y, &s, &c);
- return (CMPLX(exp_x * c, exp_x * s));
- }
+ exp_x = expl(x);
+ sincosl(y, &s, &c);
+ return (CMPLXL(exp_x * c, exp_x * s));
}
--
Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20190203173056.GA42499>