d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add cexp() and cexpf().

Browse Source 
Reviewed by:	bde (earlier version)
This commit is contained in:
das 2011-03-07 03:09:24 +00:00
parent 20bf8c09c5
commit 55e5832ebf
5 changed files with 317 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
lib/msun/Makefile
View File

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

3
lib/msun/Symbol.map
View File

@ -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
View 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
View 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
View 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)));
}
}