Add cexp() and cexpf().
Reviewed by: bde (earlier version)
This commit is contained in:
parent
799caf55bb
commit
f3732b5aaa
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=219359
@ -101,7 +101,8 @@ COMMON_SRCS+= e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
|
||||
.endif
|
||||
|
||||
# C99 complex functions
|
||||
COMMON_SRCS+= s_cimag.c s_cimagf.c s_cimagl.c s_conj.c s_conjf.c s_conjl.c \
|
||||
COMMON_SRCS+= s_cexp.c s_cexpf.c s_cimag.c s_cimagf.c s_cimagl.c \
|
||||
s_conj.c s_conjf.c s_conjl.c \
|
||||
s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_creall.c
|
||||
|
||||
# FreeBSD's C library supplies these functions:
|
||||
@ -124,7 +125,7 @@ SRCS= ${COMMON_SRCS} ${ARCH_SRCS}
|
||||
|
||||
INCS= fenv.h math.h
|
||||
|
||||
MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 ceil.3 \
|
||||
MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 ceil.3 cexp.3 \
|
||||
cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
|
||||
feclearexcept.3 feenableexcept.3 fegetenv.3 \
|
||||
fegetround.3 fenv.3 floor.3 \
|
||||
@ -143,6 +144,7 @@ MLINKS+=atanh.3 atanhf.3
|
||||
MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \
|
||||
atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3
|
||||
MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3
|
||||
MLINKS+=cexp.3 cexpf.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 \
|
||||
|
@ -222,6 +222,9 @@ FBSD_1.1 {
|
||||
/* First added in 9.0-CURRENT */
|
||||
FBSD_1.2 {
|
||||
__isnanf;
|
||||
cexp;
|
||||
cexpf;
|
||||
log2;
|
||||
log2f;
|
||||
logl;
|
||||
};
|
||||
|
113
lib/msun/man/cexp.3
Normal file
113
lib/msun/man/cexp.3
Normal file
@ -0,0 +1,113 @@
|
||||
.\" 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.
|
||||
.\"
|
||||
.\" $FreeBSD$
|
||||
.\"
|
||||
.Dd March 6, 2011
|
||||
.Dt CEXP 3
|
||||
.Os
|
||||
.Sh NAME
|
||||
.Nm cexp ,
|
||||
.Nm cexpf
|
||||
.Nd complex exponential functions
|
||||
.Sh LIBRARY
|
||||
.Lb libm
|
||||
.Sh SYNOPSIS
|
||||
.In complex.h
|
||||
.Ft double complex
|
||||
.Fn cexp "double complex z"
|
||||
.Ft float complex
|
||||
.Fn cexpf "float complex z"
|
||||
.Sh DESCRIPTION
|
||||
The
|
||||
.Fn cexp
|
||||
and
|
||||
.Fn cexpf
|
||||
functions compute the complex exponential of
|
||||
.Fa z ,
|
||||
also known as
|
||||
.Em cis Ns ( Ns
|
||||
.Fa z Ns )
|
||||
.Sh RETURN VALUES
|
||||
For real numbers
|
||||
.Fa x
|
||||
and
|
||||
.Fa y ,
|
||||
.Fn cexp
|
||||
behaves according to Euler's formula:
|
||||
.Bd -ragged -offset indent
|
||||
.Fn cexp "x + I*y"
|
||||
=
|
||||
.Ns ( Sy e Ns ** Ns
|
||||
.Fa x *
|
||||
.Em cos Ns ( Ns
|
||||
.Fa y Ns )) + ( Ns
|
||||
.Sy I
|
||||
*
|
||||
.Sy e Ns ** Ns
|
||||
.Fa x
|
||||
*
|
||||
.Em sin Ns ( Ns
|
||||
.Fa y Ns ))
|
||||
.Ed
|
||||
.Pp
|
||||
Generally speaking, infinities, zeroes and \*(Nas are handled as would
|
||||
be expected from this identity given the usual rules of floating-point
|
||||
arithmetic.
|
||||
However, care is taken to avoid generating \*(Nas when they are not deserved.
|
||||
For example, mathematically we expect that
|
||||
.Fo cimag
|
||||
.Fn cexp "x + I*0" Fc
|
||||
= 0 regardless of the value of
|
||||
.Fa x ,
|
||||
and
|
||||
.Fn cexp
|
||||
preserves this identity even if
|
||||
.Fa x
|
||||
is \*(If or \*(Na.
|
||||
Likewise,
|
||||
.Fn cexp "-\*(If + I*y"
|
||||
= 0 and
|
||||
.Fo creal
|
||||
.Fn cexp "\*(If + I*y" Fc
|
||||
= \*(If
|
||||
for any
|
||||
.Fa y
|
||||
(even though the latter property is only mathematically true for
|
||||
representable
|
||||
.Fa y . )
|
||||
If
|
||||
.Fa y
|
||||
is not finite, the sign of the result is indeterminate.
|
||||
.Sh SEE ALSO
|
||||
.Xr complex 3 ,
|
||||
.Xr exp 3 ,
|
||||
.Xr math 3 ,
|
||||
.Sh STANDARDS
|
||||
The
|
||||
.Fn cexp
|
||||
and
|
||||
.Fn cexpf
|
||||
functions conform to
|
||||
.St -isoC-99 .
|
99
lib/msun/src/s_cexp.c
Normal file
99
lib/msun/src/s_cexp.c
Normal file
@ -0,0 +1,99 @@
|
||||
/*-
|
||||
* 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"
|
||||
|
||||
static const uint32_t
|
||||
exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */
|
||||
cexp_ovfl = 0x4096b8e4, /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
|
||||
k = 1799; /* constant for reduction */
|
||||
|
||||
static const double
|
||||
kln2 = 1246.97177782734161156; /* k * ln2 */
|
||||
|
||||
double complex
|
||||
cexp(double complex z)
|
||||
{
|
||||
double x, y, exp_x;
|
||||
uint32_t hx, hy, lx, ly;
|
||||
int scale;
|
||||
|
||||
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 (cpack(exp(x), y));
|
||||
if (hy >= 0x7ff00000) {
|
||||
EXTRACT_WORDS(hx, lx, x);
|
||||
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
|
||||
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
|
||||
return (cpack(y - y, y - y));
|
||||
} else if (hx & 0x80000000) {
|
||||
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
|
||||
return (cpack(0.0, 0.0));
|
||||
} else {
|
||||
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
|
||||
return (cpack(x, y - y));
|
||||
}
|
||||
}
|
||||
|
||||
GET_HIGH_WORD(hx, x);
|
||||
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). 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 MANT_DIG to avoid loss of
|
||||
* accuracy due to underflow if sin(y) is tiny.
|
||||
*/
|
||||
exp_x = exp(x - kln2);
|
||||
GET_HIGH_WORD(hx, exp_x);
|
||||
SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 52) << 20));
|
||||
scale = (hx >> 20) - (0x3ff + 52) + k;
|
||||
return (cpack(scalbn(cos(y) * exp_x, scale),
|
||||
scalbn(sin(y) * exp_x, scale)));
|
||||
} 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);
|
||||
return (cpack(exp_x * cos(y), exp_x * sin(y)));
|
||||
}
|
||||
}
|
98
lib/msun/src/s_cexpf.c
Normal file
98
lib/msun/src/s_cexpf.c
Normal file
@ -0,0 +1,98 @@
|
||||
/*-
|
||||
* 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"
|
||||
|
||||
static const uint32_t
|
||||
exp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */
|
||||
cexp_ovfl = 0x43400074, /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */
|
||||
k = 235; /* constant for reduction */
|
||||
|
||||
static const float
|
||||
kln2 = 162.88958740f; /* k * ln2 */
|
||||
|
||||
float complex
|
||||
cexpf(float complex z)
|
||||
{
|
||||
float x, y, exp_x;
|
||||
uint32_t hx, hy;
|
||||
int scale;
|
||||
|
||||
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 (cpackf(expf(x), y));
|
||||
GET_FLOAT_WORD(hx, x);
|
||||
if (hy >= 0x7f800000) {
|
||||
if ((hx & 0x7fffffff) != 0x7f800000) {
|
||||
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
|
||||
return (cpackf(y - y, y - y));
|
||||
} else if (hx & 0x80000000) {
|
||||
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
|
||||
return (cpackf(0.0, 0.0));
|
||||
} else {
|
||||
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
|
||||
return (cpackf(x, y - y));
|
||||
}
|
||||
}
|
||||
|
||||
if (hx >= exp_ovfl && hx <= cexp_ovfl) {
|
||||
/*
|
||||
* x is between 88.7 and 192, so we must scale to avoid
|
||||
* overflow in expf(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 MANT_DIG to avoid loss of
|
||||
* accuracy due to underflow if sin(y) is tiny.
|
||||
*/
|
||||
exp_x = expf(x - kln2);
|
||||
GET_FLOAT_WORD(hx, exp_x);
|
||||
SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 23) << 23));
|
||||
scale = (hx >> 23) - (0x7f + 23) + k;
|
||||
return (cpackf(scalbnf(cosf(y) * exp_x, scale),
|
||||
scalbnf(sinf(y) * exp_x, scale)));
|
||||
} 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 = expf(x);
|
||||
return (cpackf(exp_x * cosf(y), exp_x * sinf(y)));
|
||||
}
|
||||
}
|
Loading…
Reference in New Issue
Block a user