From 55e5832ebf71c5279994dae84bf74bc646ad717f Mon Sep 17 00:00:00 2001 From: das Date: Mon, 7 Mar 2011 03:09:24 +0000 Subject: [PATCH] Add cexp() and cexpf(). Reviewed by: bde (earlier version) --- lib/msun/Makefile | 6 ++- lib/msun/Symbol.map | 3 ++ lib/msun/man/cexp.3 | 113 +++++++++++++++++++++++++++++++++++++++++ lib/msun/src/s_cexp.c | 99 ++++++++++++++++++++++++++++++++++++ lib/msun/src/s_cexpf.c | 98 +++++++++++++++++++++++++++++++++++ 5 files changed, 317 insertions(+), 2 deletions(-) create mode 100644 lib/msun/man/cexp.3 create mode 100644 lib/msun/src/s_cexp.c create mode 100644 lib/msun/src/s_cexpf.c diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 7950961d1a60..2d7d5cccb77e 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -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 \ diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index 69618e68a7d3..2ef5ebf20b2c 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -222,6 +222,9 @@ FBSD_1.1 { /* First added in 9.0-CURRENT */ FBSD_1.2 { __isnanf; + cexp; + cexpf; log2; log2f; + logl; }; diff --git a/lib/msun/man/cexp.3 b/lib/msun/man/cexp.3 new file mode 100644 index 000000000000..407284cebf84 --- /dev/null +++ b/lib/msun/man/cexp.3 @@ -0,0 +1,113 @@ +.\" Copyright (c) 2011 David Schultz +.\" 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 . diff --git a/lib/msun/src/s_cexp.c b/lib/msun/src/s_cexp.c new file mode 100644 index 000000000000..ecf0992dc21d --- /dev/null +++ b/lib/msun/src/s_cexp.c @@ -0,0 +1,99 @@ +/*- + * Copyright (c) 2011 David Schultz + * 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 +__FBSDID("$FreeBSD$"); + +#include +#include + +#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))); + } +} diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/src/s_cexpf.c new file mode 100644 index 000000000000..4ea39312c040 --- /dev/null +++ b/lib/msun/src/s_cexpf.c @@ -0,0 +1,98 @@ +/*- + * Copyright (c) 2011 David Schultz + * 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 +__FBSDID("$FreeBSD$"); + +#include +#include + +#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))); + } +}