From aaf70b231459920935d6fe1fa2e98b21770241e4 Mon Sep 17 00:00:00 2001 From: David Schultz Date: Sat, 15 Dec 2007 08:38:44 +0000 Subject: [PATCH] Implement and document csqrt(3) and csqrtf(3). --- include/complex.h | 4 +- lib/msun/Makefile | 6 ++- lib/msun/man/csqrt.3 | 94 ++++++++++++++++++++++++++++++++++++ lib/msun/src/s_csqrt.c | 104 ++++++++++++++++++++++++++++++++++++++++ lib/msun/src/s_csqrtf.c | 88 ++++++++++++++++++++++++++++++++++ 5 files changed, 293 insertions(+), 3 deletions(-) create mode 100644 lib/msun/man/csqrt.3 create mode 100644 lib/msun/src/s_csqrt.c create mode 100644 lib/msun/src/s_csqrtf.c diff --git a/include/complex.h b/include/complex.h index 5627f61f2653..24251f20f7ab 100644 --- a/include/complex.h +++ b/include/complex.h @@ -1,5 +1,5 @@ /*- - * Copyright (c) 2001 The FreeBSD Project. + * Copyright (c) 2001-2007 The FreeBSD Project. * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -57,6 +57,8 @@ long double complex double creal(double complex); float crealf(float complex); long double creall(long double complex); +double complex csqrt(double complex); +float complex csqrtf(float complex); __END_DECLS diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 8a04ef209035..eaad0f6c3d0c 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -40,7 +40,8 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \ k_tan.c k_tanf.c \ s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c \ s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c s_ceill.c \ - s_copysign.c s_copysignf.c s_cos.c s_cosf.c s_erf.c s_erff.c \ + s_copysign.c s_copysignf.c s_cos.c s_cosf.c \ + s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \ s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \ s_finite.c s_finitef.c \ s_floor.c s_floorf.c s_floorl.c s_fma.c s_fmaf.c \ @@ -93,7 +94,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 \ - cimag.3 copysign.3 cos.3 cosh.3 erf.3 exp.3 fabs.3 fdim.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 \ fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \ @@ -114,6 +115,7 @@ MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \ MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3 MLINKS+=cos.3 cosf.3 MLINKS+=cosh.3 coshf.3 +MLINKS+=csqrt.3 csqrtf.3 MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 MLINKS+=exp.3 expm1.3 exp.3 log.3 exp.3 log10.3 exp.3 log1p.3 exp.3 pow.3 \ exp.3 exp2.3 exp.3 exp2f.3 exp.3 expf.3 \ diff --git a/lib/msun/man/csqrt.3 b/lib/msun/man/csqrt.3 new file mode 100644 index 000000000000..132b314394b9 --- /dev/null +++ b/lib/msun/man/csqrt.3 @@ -0,0 +1,94 @@ +.\" Copyright (c) 2007 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 December 14, 2007 +.Dt CSQRT 3 +.Os +.Sh NAME +.Nm csqrt , +.Nm csqrtf +.Nd complex square root functions +.Sh LIBRARY +.Lb libm +.Sh SYNOPSIS +.In complex.h +.Ft double complex +.Fn csqrt "double complex z" +.Ft float complex +.Fn csqrtf "float complex z" +.Sh DESCRIPTION +The +.Fn csqrt +and +.Fn csqrtf +functions compute the square root of +.Fa z +in the complex plane, with a branch cut along the negative real axis. +In other words, +.Fn csqrt +and +.Fn csqrtf +always return the square root whose real part is non-negative. +.Sh RETURN VALUES +These functions return the requested square root. +The square root of 0 is +.Li +0 \*(Pm 0 , +where the imaginary parts of the input and respective result have +the same sign. +For infinities and \*(Nas, the following rules apply, with the +earlier rules having precedence: +.Bl -column -offset indent "-\*(If + \*(Na*I" "\*(If \*(Pm \*(If*I " "(for all k)" +.Em Input Result +k + \*(If*I \*(If + \*(If*I (for all k) +-\*(If + \*(Na*I \*(Na \*(Pm \*(If*I +\*(If + \*(Na*I \*(If + \*(Na*I +k + \*(Na*I \*(Na + \*(Na*I +\*(Na + k*I \*(Na + \*(Na*I +-\*(If + k*I +0 + \*(If*I +\*(If + k*I \*(If + 0*I +.El +.Pp +For numbers with negative imaginary parts, the above special cases +apply given the identity: +.Dl csqrt(conj(z) = conj(sqrt(z)) +Note that the sign of \*(Na is indeterminate. +Also, if the real or imaginary part of the input is finite and +an \*(Na is generated, an invalid exception will be thrown. +.Sh SEE ALSO +.Xr cabs 3 , +.Xr fenv 3 , +.Xr math 3 , +.Sh STANDARDS +The +.Fn csqrt +and +.Fn csqrtf +functions conform to +.St -isoC-99 . +.Sh BUGS +For +.Fn csqrt , +inexact results are not correctly rounded in general. diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c new file mode 100644 index 000000000000..3b1e18721938 --- /dev/null +++ b/lib/msun/src/s_csqrt.c @@ -0,0 +1,104 @@ +/*- + * Copyright (c) 2007 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" + +/* + * gcc doesn't implement complex multiplication or division correctly, + * so we need to handle infinities specially. We turn on this pragma to + * notify conforming c99 compilers that the fast-but-incorrect code that + * gcc generates is acceptable, since the special cases have already been + * handled. + */ +#pragma STDC CX_LIMITED_RANGE on + +/* We risk spurious overflow for components >= DBL_MAX/(1+sqrt(2)) */ +#define THRESH 0x1.a827999fcef32p+1022 + +double complex +csqrt(double complex z) +{ + double a = creal(z), b = cimag(z); + double t; + double complex result; + int scale; + + /* Handle special cases. */ + if (z == 0) + return (cpack(0, b)); + if (isinf(b)) + return (cpack(INFINITY, b)); + if (isnan(a)) { + t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ + return (cpack(t, t)); /* return NaN + NaN i */ + } + if (isinf(a)) { + /* + * csqrt(inf + nan i) = inf + nan i + * csqrt(inf + y i) = inf + 0 i + * csqrt(-inf + nan i) = nan +- inf i + * csqrt(-inf + y i) = 0 + inf i + */ + if (signbit(a)) + return (cpack(fabs(b - b), copysign(a, b))); + else + return (cpack(a, copysign(b - b, b))); + } + /* + * The remaining special case (b is NaN) is handled just fine by + * the normal code path below. + */ + + /* Scale to avoid overflow. */ + if (a >= THRESH || b >= THRESH) { + a *= 0.25; + b *= 0.25; + scale = 1; + } else { + scale = 0; + } + + /* Algorithm 312, CACM vol 10, Oct 1967 */ + if (a >= 0) { + t = sqrt((a + hypot(a, b)) * 0.5); + result = cpack(t, b / (2 * t)); + } else { + t = sqrt((-a + hypot(a, b)) * 0.5); + result = cpack(fabs(b) / (2 * t), copysign(t, b)); + } + + /* Rescale */ + if (scale) + return (result * 2); + else + return (result); +} diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c new file mode 100644 index 000000000000..cc5485a5b23a --- /dev/null +++ b/lib/msun/src/s_csqrtf.c @@ -0,0 +1,88 @@ +/*- + * Copyright (c) 2007 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" + +/* + * gcc doesn't implement complex multiplication or division correctly, + * so we need to handle infinities specially. We turn on this pragma to + * notify conforming c99 compilers that the fast-but-incorrect code that + * gcc generates is acceptable, since the special cases have already been + * handled. + */ +#pragma STDC CX_LIMITED_RANGE on + +float complex +csqrtf(float complex z) +{ + float a = crealf(z), b = cimagf(z); + double t; + + /* Handle special cases. */ + if (z == 0) + return (cpackf(0, b)); + if (isinf(b)) + return (cpackf(INFINITY, b)); + if (isnan(a)) { + t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ + return (cpackf(t, t)); /* return NaN + NaN i */ + } + if (isinf(a)) { + /* + * csqrtf(inf + nan i) = inf + nan i + * csqrtf(inf + y i) = inf + 0 i + * csqrtf(-inf + nan i) = nan +- inf i + * csqrtf(-inf + y i) = 0 + inf i + */ + if (signbit(a)) + return (cpackf(fabsf(b - b), copysignf(a, b))); + else + return (cpackf(a, copysignf(b - b, b))); + } + /* + * The remaining special case (b is NaN) is handled just fine by + * the normal code path below. + */ + + /* + * We compute t in double precision to avoid overflow and to + * provide correct rounding in nearly all cases. + * This is Algorithm 312, CACM vol 10, Oct 1967. + */ + if (a >= 0) { + t = sqrt((a + hypot(a, b)) * 0.5); + return (cpackf(t, b / (2.0 * t))); + } else { + t = sqrt((-a + hypot(a, b)) * 0.5); + return (cpackf(fabsf(b) / (2.0 * t), copysignf(t, b))); + } +}