Attachment #202398 for bug #235413

Lines 243-248 Link Here

(-)lib/msun/src/math_private.h (+18 lines)
243	} while (0)	243	} while (0)
244		244
245	/*	245	/*
		246	* Get expsign and mantissa as 16 bit and 2 32-bit ints from an 80 bit long
		247	* double.
		248	*/
		249
		250	#define EXTRACT_LDBL80_WORDS2(ix0,ix1,ix2,d) \
		251	do { \
		252	union IEEEl2bits ew_u; \
		253	ew_u.e = (d); \
		254	(ix0) = ew_u.xbits.expsign; \
		255	(ix1) = ew_u.bits.manh; \
		256	(ix2) = ew_u.bits.manl; \
		257	} while (0)
		258
		259	/*
246	* Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit	260	* Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
247	* long double.	261	* long double.
248	*/	262	*/
Lines 920-924 Link Here
920	long double __kernel_sinl(long double, long double, int);	934	long double __kernel_sinl(long double, long double, int);
921	long double __kernel_cosl(long double, long double);	935	long double __kernel_cosl(long double, long double);
922	long double __kernel_tanl(long double, long double, int);	936	long double __kernel_tanl(long double, long double, int);
		937	long double __ldexp_expl(long double, int);
		938	#ifdef _COMPLEX_H
		939	long double complex __ldexp_cexpl(long double complex, int);
		940	#endif
923		941
924	#endif /* !_MATH_PRIVATE_H_ */	942	#endif /* !_MATH_PRIVATE_H_ */




__FBSDID("$FreeBSD$");

#include <complex.h>
#include <float.h>
#include <math.h>

#include "math_private.h"

double complex
cexp(double complex z)
{
	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;


	/* 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) {

		 *  -  x = NaN (spurious inexact exception from y)
		 */
		exp_x = exp(x);
		sincos(y, &s, &c);
		return (CMPLX(exp_x * c, exp_x * s));
	}
}

#if (LDBL_MANT_DIG == 53)
__weak_reference(cexp, cexpl);
#endif




float complex
cexpf(float complex z)
{
	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;


	/* 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) {

		 *  -  x = NaN (spurious inexact exception from y)
		 */
		exp_x = expf(x);
		sincosf(y, &s, &c);
		return (CMPLXF(exp_x * c, exp_x * s));
	}
}

Lines 88-94 Link Here

(-)lib/msun/src/k_exp.c (-3 / +4 lines)
88	double complex	88	double complex
89	__ldexp_cexp(double complex z, int expt)	89	__ldexp_cexp(double complex z, int expt)
90	{	90	{
91	double x, y, exp_x, scale1, scale2;	91	double c, exp_x, s, scale1, scale2, x, y;
92	int ex_expt, half_expt;	92	int ex_expt, half_expt;
93		93
94	x = creal(z);	94	x = creal(z);
Lines 105-110 Link Here
105	half_expt = expt - half_expt;	105	half_expt = expt - half_expt;
106	INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);	106	INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
107		107
108	return (CMPLX(cos(y) * exp_x * scale1 * scale2,	108	sincos(y, &s, &c);
109	sin(y) * exp_x * scale1 * scale2));	109	return (CMPLX(c * exp_x * scale1 * scale2,
		110	s * exp_x * scale1 * scale2));
110	}	111	}

Lines 71-77 Link Here

(-)lib/msun/src/k_expf.c (-3 / +4 lines)
71	float complex	71	float complex
72	__ldexp_cexpf(float complex z, int expt)	72	__ldexp_cexpf(float complex z, int expt)
73	{	73	{
74	float x, y, exp_x, scale1, scale2;	74	float c, exp_x, s, scale1, scale2, x, y;
75	int ex_expt, half_expt;	75	int ex_expt, half_expt;
76		76
77	x = crealf(z);	77	x = crealf(z);
Lines 84-89 Link Here
84	half_expt = expt - half_expt;	84	half_expt = expt - half_expt;
85	SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);	85	SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
86		86
87	return (CMPLXF(cosf(y) * exp_x * scale1 * scale2,	87	sincosf(y, &s, &c);
88	sinf(y) * exp_x * scale1 * scale2));	88	return (CMPLXF(c * exp_x * scale1 * scale2,
		89	s * exp_x * scale1 * scale2));
89	}	90	}




/*-
 * 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);
}




/*-
 * 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");




/*-
 * 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));
}




/*-
 * 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));
	}
}




.\"
.\" $FreeBSD$
.\"
.Dd February 4, 2019
.Dt CEXP 3
.Os
.Sh NAME
.Nm cexp ,
.Nm cexpf ,
.Nm cexpl
.Nd complex exponential functions
.Sh LIBRARY
.Lb libm

.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 cexpf ,
and
.Fn cexpl
functions compute the complex exponential of
.Fa z ,
also known as

.Xr math 3
.Sh STANDARDS
The
.Fn cexp ,
.Fn cexpf ,
and
.Fn cexpl
functions conform to
.St -isoC-99 .

Lines 98-103 Link Here

(-)include/complex.h (+2 lines)
98	float complex ccoshf(float complex);	98	float complex ccoshf(float complex);
99	double complex cexp(double complex);	99	double complex cexp(double complex);
100	float complex cexpf(float complex);	100	float complex cexpf(float complex);
		101	long double complex
		102	cexpl(long double complex);
101	double cimag(double complex) __pure2;	103	double cimag(double complex) __pure2;
102	float cimagf(float complex) __pure2;	104	float cimagf(float complex) __pure2;
103	long double cimagl(long double complex) __pure2;	105	long double cimagl(long double complex) __pure2;

Return to bug 235413

Lines 30-35 Link Here

(-)lib/msun/src/s_cexp.c (-6 / +15 lines)
30	__FBSDID("$FreeBSD$");	30	__FBSDID("$FreeBSD$");
31		31
32	#include <complex.h>	32	#include <complex.h>
		33	#include <float.h>
33	#include <math.h>	34	#include <math.h>
34		35
35	#include "math_private.h"	36	#include "math_private.h"
Lines 41-52 Link Here
41	double complex	42	double complex
42	cexp(double complex z)	43	cexp(double complex z)
43	{	44	{
44	double x, y, exp_x;	45	double c, exp_x, s, x, y;
45	uint32_t hx, hy, lx, ly;	46	uint32_t hx, hy, lx, ly;
46		47
47	x = creal(z);	48	x = creal(z);
48	y = cimag(z);	49	y = cimag(z);
49		50
		51	EXTRACT_WORDS(hx, lx, x);
		52	/* cexp(0 + I y) = cos(y) + I sin(y) */
		53	if (((hx & 0x7fffffff) \| lx) == 0) {
		54	sincos(y, &s, &c);
		55	return (CMPLX(c, s));
		56	}
		57
50	EXTRACT_WORDS(hy, ly, y);	58	EXTRACT_WORDS(hy, ly, y);
51	hy &= 0x7fffffff;	59	hy &= 0x7fffffff;
52		60
Lines 53-62 Link Here
53	/* cexp(x + I 0) = exp(x) + I 0 */	61	/* cexp(x + I 0) = exp(x) + I 0 */
54	if ((hy \| ly) == 0)	62	if ((hy \| ly) == 0)
55	return (CMPLX(exp(x), y));	63	return (CMPLX(exp(x), y));
56	EXTRACT_WORDS(hx, lx, x);
57	/* cexp(0 + I y) = cos(y) + I sin(y) */
58	if (((hx & 0x7fffffff) \| lx) == 0)
59	return (CMPLX(cos(y), sin(y)));
60		64
61	if (hy >= 0x7ff00000) {	65	if (hy >= 0x7ff00000) {
62	if (lx != 0 \|\| (hx & 0x7fffffff) != 0x7ff00000) {	66	if (lx != 0 \|\| (hx & 0x7fffffff) != 0x7ff00000) {
Lines 86-91 Link Here
86	* - x = NaN (spurious inexact exception from y)	90	* - x = NaN (spurious inexact exception from y)
87	*/	91	*/
88	exp_x = exp(x);	92	exp_x = exp(x);
89	return (CMPLX(exp_x * cos(y), exp_x * sin(y)));	93	sincos(y, &s, &c);
		94	return (CMPLX(exp_x * c, exp_x * s));
90	}	95	}
91	}	96	}
		97
		98	#if (LDBL_MANT_DIG == 53)
		99	__weak_reference(cexp, cexpl);
		100	#endif

Lines 41-52 Link Here

(-)lib/msun/src/s_cexpf.c (-6 / +10 lines)
41	float complex	41	float complex
42	cexpf(float complex z)	42	cexpf(float complex z)
43	{	43	{
44	float x, y, exp_x;	44	float c, exp_x, s, x, y;
45	uint32_t hx, hy;	45	uint32_t hx, hy;
46		46
47	x = crealf(z);	47	x = crealf(z);
48	y = cimagf(z);	48	y = cimagf(z);
49		49
		50	GET_FLOAT_WORD(hx, x);
		51	/* cexp(0 + I y) = cos(y) + I sin(y) */
		52	if ((hx & 0x7fffffff) == 0) {
		53	sincosf(y, &s, &c);
		54	return (CMPLXF(c, s));
		55	}
		56
50	GET_FLOAT_WORD(hy, y);	57	GET_FLOAT_WORD(hy, y);
51	hy &= 0x7fffffff;	58	hy &= 0x7fffffff;
52		59
Lines 53-62 Link Here
53	/* cexp(x + I 0) = exp(x) + I 0 */	60	/* cexp(x + I 0) = exp(x) + I 0 */
54	if (hy == 0)	61	if (hy == 0)
55	return (CMPLXF(expf(x), y));	62	return (CMPLXF(expf(x), y));
56	GET_FLOAT_WORD(hx, x);
57	/* cexp(0 + I y) = cos(y) + I sin(y) */
58	if ((hx & 0x7fffffff) == 0)
59	return (CMPLXF(cosf(y), sinf(y)));
60		63
61	if (hy >= 0x7f800000) {	64	if (hy >= 0x7f800000) {
62	if ((hx & 0x7fffffff) != 0x7f800000) {	65	if ((hx & 0x7fffffff) != 0x7f800000) {
Lines 86-91 Link Here
86	* - x = NaN (spurious inexact exception from y)	89	* - x = NaN (spurious inexact exception from y)
87	*/	90	*/
88	exp_x = expf(x);	91	exp_x = expf(x);
89	return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));	92	sincosf(y, &s, &c);
		93	return (CMPLXF(exp_x * c, exp_x * s));
90	}	94	}
91	}	95	}

Line 0 Link Here

(-)lib/msun/ld128/k_cexpl.c (+126 lines)
		1	/*-
		2	* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
		3	*
		4	* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
		5	* All rights reserved.
		6	*
		7	* Redistribution and use in source and binary forms, with or without
		8	* modification, are permitted provided that the following conditions
		9	* are met:
		10	* 1. Redistributions of source code must retain the above copyright
		11	* notice, this list of conditions and the following disclaimer.
		12	* 2. Redistributions in binary form must reproduce the above copyright
		13	* notice, this list of conditions and the following disclaimer in the
		14	* documentation and/or other materials provided with the distribution.
		15	*
		16	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
		17	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		18	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		19	* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
		20	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		21	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		22	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		23	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		24	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		25	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		26	* SUCH DAMAGE.
		27	*/
		28
		29	#include <sys/cdefs.h>
		30	__FBSDID("$FreeBSD$");
		31
		32	#include <complex.h>
		33
		34	#include "math.h"
		35	#include "math_private.h"
		36
		37	#warning These functions are broken on ld128 architectures.
		38	#warning Functions defined here are currently unused.
		39	#warning Someone who cares needs to convert src/k_exp.c.
		40
		41	#if 0
		42	static const uint32_t k = 1799; /* constant for reduction */
		43	static const double kln2 = 1246.97177782734161156; /* k * ln2 */
		44	#endif
		45
		46	/*
		47	* Compute exp(x), scaled to avoid spurious overflow. An exponent is
		48	* returned separately in 'expt'.
		49	*
		50	* Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
		51	* Output: 21023 <= y < 21024
		52	*/
		53	static long double
		54	__frexp_expl(long double x, int *expt)
		55	{
		56	#if 0
		57	double exp_x;
		58	uint32_t hx;
		59
		60	/*
		61	* We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
		62	* minimize \|exp(kln2) - 2**k\|. We also scale the exponent of
		63	* exp_x to MAX_EXP so that the result can be multiplied by
		64	* a tiny number without losing accuracy due to denormalization.
65	*/
66	exp_x = exp(x - kln2);
67	GET_HIGH_WORD(hx, exp_x);
68	*expt = (hx >> 20) - (0x3ff + 1023) + k;
69	SET_HIGH_WORD(exp_x, (hx & 0xfffff) \| ((0x3ff + 1023) << 20));
70	return (exp_x);
71	#endif
72	return (x - x) / (x - x);
73	}
74
75	/*
76	* __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
77	* They are intended for large arguments (real part >= ln(DBL_MAX))
78	* where care is needed to avoid overflow.
79	*
80	* The present implementation is narrowly tailored for our hyperbolic and
81	* exponential functions. We assume expt is small (0 or -1), and the caller
82	* has filtered out very large x, for which overflow would be inevitable.
83	*/
84
85	long double
86	__ldexp_expl(long double x, int expt)
87	{
88	#if 0
89	double exp_x, scale;
90	int ex_expt;
91
92	exp_x = __frexp_exp(x, &ex_expt);
93	expt += ex_expt;
94	INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);
95	return (exp_x * scale);
96	#endif
97	return (x - x) / (x - x);
98	}
99
100	long double complex
101	__ldexp_cexpl(long double complex z, int expt)
102	{
103	#if 0
104	double x, y, exp_x, scale1, scale2;
105	int ex_expt, half_expt;
106
107	x = creal(z);
108	y = cimag(z);
109	exp_x = __frexp_exp(x, &ex_expt);
110	expt += ex_expt;
111
112	/*
113	* Arrange so that scale1 * scale2 == 2**expt. We use this to
114	* compensate for scalbn being horrendously slow.
115	*/
116	half_expt = expt / 2;
117	INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0);
118	half_expt = expt - half_expt;
119	INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
120
121	sincos(y, &s, &c);
122	return (CMPLX(c) * exp_x * scale1 * scale2,
123	s * exp_x * scale1 * scale2));
124	#endif
125	return (x - x) / (x - x);
126	}

Line 0 Link Here

(-)lib/msun/ld128/s_cexpl.c (+76 lines)
		1	/*-
		2	* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
		3	*
		4	* Copyright (c) 2019 Steven G. Kargl <kargl@FreeBSD.ORG>
		5	* All rights reserved.
		6	*
		7	* Redistribution and use in source and binary forms, with or without
		8	* modification, are permitted provided that the following conditions
		9	* are met:
		10	* 1. Redistributions of source code must retain the above copyright
		11	* notice, this list of conditions and the following disclaimer.
		12	* 2. Redistributions in binary form must reproduce the above copyright
		13	* notice, this list of conditions and the following disclaimer in the
		14	* documentation and/or other materials provided with the distribution.
		15	*
		16	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
		17	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		18	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		19	* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
		20	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		21	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		22	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		23	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		24	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		25	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		26	* SUCH DAMAGE.
		27	*/
		28
		29	#include <sys/cdefs.h>
		30	__FBSDID("$FreeBSD$");
		31
		32	#include <complex.h>
		33	#include <float.h>
		34	#include <math.h>
		35
		36	#include "math_private.h"
		37
		38	#warning cexpl is likely broken on ld128 architectures.
		39	#warning Someone who cares needs to convert src/s_cexp.c.
		40
		41	long double complex
		42	cexpl(long double complex z)
		43	{
		44	long double c, exp_x, s, x, y;
		45
		46	x = creall(z);
		47	y = cimagl(z);
		48
		49	/* cexp(0 + I y) = cos(y) + I sin(y) */
		50	if (x == 0) {
		51	sincosl(y, &s, &c);
		52	return (CMPLXL(c, s));
		53	}
		54
		55	/* cexp(x + I 0) = exp(x) + I 0 */
		56	if (y == 0)
		57	return (CMPLXL(expl(x), y));
		58
		59	if (!isfinite(y)) {
		60	if (isfinite(x) \|\| isnan(x)) {
		61	/* cexp(finite\|NaN +- I Inf\|NaN) = NaN + I NaN */
		62	return (CMPLXL(y - y, y - y));
		63	} else if (isinf(x) && copysignl(1.L, x) < 0) {
		64	/* cexp(-Inf +- I Inf\|NaN) = 0 + I 0 */
65	return (CMPLXL(0.0L, 0.0L));
66	} else {
67	/* cexp(+Inf +- I Inf\|NaN) = Inf + I NaN */
68	return (CMPLXL(x, y - y));
69	}
70	}
71
72	exp_x = expl(x);
73	sincosl(y, &s, &c);
74	return (CMPLXL(exp_x * c, exp_x * s));
75	}
76	__warn_references("Precision of cexpl may have lower than expected precision");

Line 0 Link Here

(-)lib/msun/ld80/k_cexpl.c (+118 lines)
		1	/*-
		2	* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
		3	*
		4	* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
		5	* All rights reserved.
		6	*
		7	* Redistribution and use in source and binary forms, with or without
		8	* modification, are permitted provided that the following conditions
		9	* are met:
		10	* 1. Redistributions of source code must retain the above copyright
		11	* notice, this list of conditions and the following disclaimer.
		12	* 2. Redistributions in binary form must reproduce the above copyright
		13	* notice, this list of conditions and the following disclaimer in the
		14	* documentation and/or other materials provided with the distribution.
		15	*
		16	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
		17	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		18	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		19	* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
		20	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		21	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		22	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		23	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		24	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		25	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		26	* SUCH DAMAGE.
		27	*
		28	* src/k_exp.c converted to long double complex by Steven G. Kargl
		29	*/
		30
		31	#include <sys/cdefs.h>
		32	__FBSDID("$FreeBSD$");
		33
		34	#include <complex.h>
		35
		36	#include "fpmath.h"
		37	#include "math.h"
		38	#include "math_private.h"
		39
		40	static const uint32_t k = 22956; /* constant for reduction */
		41	static const union IEEEl2bits
		42	kln2u = LD80C(0xf89f8bf509c862ca, 13, 1.59118866769341045231e+04L);
		43	#define kln2 (kln2u.e)
		44	/*
		45	* Compute exp(x), scaled to avoid spurious overflow. An exponent is
		46	* returned separately in 'expt'.
		47	*
		48	* Input: ln(LDBL_MAX) <= x < ln(2 * LDBL_MAX / LDBL_MIN_DENORM) ~= 22756.0
		49	* Output: 216383 <= y < 216384
		50	*/
		51	static long double
		52	__frexp_expl(long double x, int *expt)
		53	{
		54	union IEEEl2bits exp_x;
		55
		56	/*
		57	* We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
		58	* minimize \|exp(kln2) - 2**k\|. We also scale the exponent of
		59	* exp_x to MAX_EXP so that the result can be multiplied by
		60	* a tiny number without losing accuracy due to denormalization.
		61	*/
		62	exp_x.e = expl(x - kln2);
		63	*expt = exp_x.bits.exp - (0x3fff + 16383) + k - 1;
		64	exp_x.bits.exp = 0x3fff + 16383;
65	return (exp_x.e);
66	}
67
68	/*
69	* __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
70	* They are intended for large arguments (real part >= ln(DBL_MAX))
71	* where care is needed to avoid overflow.
72	*
73	* The present implementation is narrowly tailored for our hyperbolic and
74	* exponential functions. We assume expt is small (0 or -1), and the caller
75	* has filtered out very large x, for which overflow would be inevitable.
76	*/
77
78	long double
79	__ldexp_expl(long double x, int expt)
80	{
81	union IEEEl2bits scale;
82	long double exp_x;
83	int ex_expt;
84
85	exp_x = __frexp_expl(x, &ex_expt);
86	expt += ex_expt;
87	scale.e = 1;
88	scale.bits.exp = 0x3fff + expt;
89	return (exp_x * scale.e);
90	}
91
92	long double complex
93	__ldexp_cexpl(long double complex z, int expt)
94	{
95	union IEEEl2bits scale1, scale2;
96	long double c, exp_x, s, x, y;
97	int ex_expt, half_expt;
98
99	x = creall(z);
100	y = cimagl(z);
101	exp_x = __frexp_expl(x, &ex_expt);
102	expt += ex_expt;
103
104	/*
105	* Arrange so that scale1 * scale2 == 2**expt. We use this to
106	* compensate for scalbn being horrendously slow.
107	*/
108	half_expt = expt / 2;
109	scale1.e = 1;
110	scale1.bits.exp = 0x3fff + half_expt;
111	half_expt = expt - half_expt;
112	scale2.e = 1;
113	scale2.bits.exp = 0x3fff + half_expt + 1;
114
115	sincosl(y, &s, &c);
116	return (CMPLXL(c * exp_x * scale1.e * scale2.e,
117	s * exp_x * scale1.e * scale2.e));
118	}

Line 0 Link Here

(-)lib/msun/ld80/s_cexpl.c (+105 lines)
		1	/*-
		2	* SPDX-License-Identifier: BSD-2-Clause-FreeBSD
		3	*
		4	* Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
		5	* All rights reserved.
		6	*
		7	* Redistribution and use in source and binary forms, with or without
		8	* modification, are permitted provided that the following conditions
		9	* are met:
		10	* 1. Redistributions of source code must retain the above copyright
		11	* notice, this list of conditions and the following disclaimer.
		12	* 2. Redistributions in binary form must reproduce the above copyright
		13	* notice, this list of conditions and the following disclaimer in the
		14	* documentation and/or other materials provided with the distribution.
		15	*
		16	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
		17	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		18	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		19	* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
		20	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		21	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		22	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		23	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		24	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		25	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		26	* SUCH DAMAGE.
		27	*
		28	* src/s_cexp.c converted to long double complex by Steven G. Kargl
		29	*/
		30
		31	#include <sys/cdefs.h>
		32	__FBSDID("$FreeBSD$");
		33
		34	#include <complex.h>
		35	#include <float.h>
		36	#ifdef __i386__
		37	#include <ieeefp.h>
		38	#endif
		39	#include <math.h>
		40
		41	#include "fpmath.h"
		42	#include "math_private.h"
		43
		44	static const uint32_t
		45	exp_ovfl = 0xb17217f7, /* high bits of MAX_EXP * ln2 ~= 11356 */
		46	cexp_ovfl = 0xb1c6a857; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
		47
		48	long double complex
		49	cexpl (long double complex z)
		50	{
		51	long double c, exp_x, s, x, y;
		52	uint32_t hwx, hwy, lwx, lwy;
		53	uint16_t hx, hy;
		54
		55	ENTERI(z);
		56
		57	x = creall(z);
		58	y = cimagl(z);
		59
		60	EXTRACT_LDBL80_WORDS2(hx, hwx, lwx, x);
		61	EXTRACT_LDBL80_WORDS2(hy, hwy, lwy, y);
		62
		63	/* cexp(0 + I y) = cos(y) + I sin(y) */
		64	if ((hwx \| lwx \| (hx & 0x7fff)) == 0) {
65	sincosl(y, &s, &c);
66	RETURNI(CMPLXL(c, s));
67	}
68
69	/* cexp(x + I 0) = exp(x) + I 0 */
70	if ((hwy \| lwy \| (hy & 0x7fff)) == 0)
71	RETURNI(CMPLXL(expl(x), y));
72
73	if (hy >= 0x7fff) {
74	if (hx < 0x7fff \|\|
75	(hx == 0x7fff && (hwx & 0x7fffffff) != 0)) {
76	/* cexp(finite\|NaN +- I Inf\|NaN) = NaN + I NaN */
77	RETURNI(CMPLXL(y - y, y - y));
78	} else if (hx & 0x8000) {
79	/* cexp(-Inf +- I Inf\|NaN) = 0 + I 0 */
80	RETURNI(CMPLXL(0.0L, 0.0L));
81	} else {
82	/* cexp(+Inf +- I Inf\|NaN) = Inf + I NaN */
83	RETURNI(CMPLXL(x, y - y));
84	}
85	}
86
87	if (hx >= 0x400c && hx <= 0x400d) {
88	/*
89	* x is between 11356 and 22755, so we must scale to avoid
90	* overflow in exp(x).
91	*/
92	RETURNI(__ldexp_cexpl(z, 0));
93	} else {
94	/*
95	* Cases covered here:
96	* - x < exp_ovfl and exp(x) won't overflow (common case)
97	* - x > cexp_ovfl, so exp(x) * s overflows for all s > 0
98	* - x = +-Inf (generated by exp())
99	* - x = NaN (spurious inexact exception from y)
100	*/
101	exp_x = expl(x);
102	sincosl(y, &s, &c);
103	RETURNI(CMPLXL(exp_x * c, exp_x * s));
104	}
105	}

Lines 24-35 Link Here

(-)lib/msun/man/cexp.3 (-6 / +11 lines)
24	.\"	24	.\"
25	.\" $FreeBSD$	25	.\" $FreeBSD$
26	.\"	26	.\"
27	.Dd March 6, 2011	27	.Dd February 4, 2019
28	.Dt CEXP 3	28	.Dt CEXP 3
29	.Os	29	.Os
30	.Sh NAME	30	.Sh NAME
31	.Nm cexp ,	31	.Nm cexp ,
32	.Nm cexpf	32	.Nm cexpf ,
		33	.Nm cexpl
33	.Nd complex exponential functions	34	.Nd complex exponential functions
34	.Sh LIBRARY	35	.Sh LIBRARY
35	.Lb libm	36	.Lb libm
Lines 39-49 Link Here
39	.Fn cexp "double complex z"	40	.Fn cexp "double complex z"
40	.Ft float complex	41	.Ft float complex
41	.Fn cexpf "float complex z"	42	.Fn cexpf "float complex z"
		43	.Ft long double complex
		44	.Fn cexpl "long double complex z"
42	.Sh DESCRIPTION	45	.Sh DESCRIPTION
43	The	46	The
44	.Fn cexp	47	.Fn cexp ,
		48	.Fn cexpf ,
45	and	49	and
46	.Fn cexpf	50	.Fn cexpl
47	functions compute the complex exponential of	51	functions compute the complex exponential of
48	.Fa z ,	52	.Fa z ,
49	also known as	53	also known as
Lines 106-113 Link Here
106	.Xr math 3	110	.Xr math 3
107	.Sh STANDARDS	111	.Sh STANDARDS
108	The	112	The
109	.Fn cexp	113	.Fn cexp ,
		114	.Fn cexpf ,
110	and	115	and
111	.Fn cexpf	116	.Fn cexpl
112	functions conform to	117	functions conform to
113	.St -isoC-99 .	118	.St -isoC-99 .