From 7e1a7394d99ee57106e7edda77a38a6cd01ecf9c Mon Sep 17 00:00:00 2001 From: das Date: Sun, 30 Mar 2008 20:07:15 +0000 Subject: Implement csqrtl(). --- lib/msun/Makefile | 4 +- lib/msun/Symbol.map | 1 + lib/msun/man/csqrt.3 | 30 +++++++++----- lib/msun/src/s_csqrt.c | 5 +++ lib/msun/src/s_csqrtl.c | 108 ++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 135 insertions(+), 13 deletions(-) create mode 100644 lib/msun/src/s_csqrtl.c (limited to 'lib/msun') diff --git a/lib/msun/Makefile b/lib/msun/Makefile index c208261..9b2e8fe 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -80,7 +80,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c .if ${LDBL_PREC} != 53 # If long double != double use these; otherwise, we alias the double versions. COMMON_SRCS+= e_hypotl.c e_sqrtl.c k_cosl.c k_sinl.c k_tanl.c \ - s_ceill.c s_cosl.c s_exp2l.c s_floorl.c s_fmal.c \ + s_ceill.c s_cosl.c s_csqrtl.c s_exp2l.c s_floorl.c s_fmal.c \ s_frexpl.c s_logbl.c s_nanl.c s_nextafterl.c s_nexttoward.c \ s_rintl.c s_scalbnl.c s_sinl.c s_tanl.c s_truncl.c w_cabsl.c .endif @@ -127,7 +127,7 @@ MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \ MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3 MLINKS+=cos.3 cosf.3 cos.3 cosl.3 MLINKS+=cosh.3 coshf.3 -MLINKS+=csqrt.3 csqrtf.3 +MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3 MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 pow.3 exp.3 powf.3 \ exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3 exp.3 expf.3 diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index e0fe88a..42258f1 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -205,4 +205,5 @@ FBSD_1.1 { sqrtl; hypotl; cabsl; + csqrtl; }; diff --git a/lib/msun/man/csqrt.3 b/lib/msun/man/csqrt.3 index 132b314..ef5445f 100644 --- a/lib/msun/man/csqrt.3 +++ b/lib/msun/man/csqrt.3 @@ -1,4 +1,4 @@ -.\" Copyright (c) 2007 David Schultz +.\" Copyright (c) 2007-2008 David Schultz .\" All rights reserved. .\" .\" Redistribution and use in source and binary forms, with or without @@ -24,12 +24,13 @@ .\" .\" $FreeBSD$ .\" -.Dd December 14, 2007 +.Dd March 30, 2008 .Dt CSQRT 3 .Os .Sh NAME .Nm csqrt , -.Nm csqrtf +.Nm csqrtf , +.Nm csqrtl .Nd complex square root functions .Sh LIBRARY .Lb libm @@ -39,18 +40,22 @@ .Fn csqrt "double complex z" .Ft float complex .Fn csqrtf "float complex z" +.Ft long double complex +.Fn csqrtl "long double complex z" .Sh DESCRIPTION The -.Fn csqrt +.Fn csqrt , +.Fn csqrtf , and -.Fn csqrtf +.Fn csqrtl 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 +.Fn csqrt , +.Fn csqrtf , and -.Fn csqrtf +.Fn csqrtl always return the square root whose real part is non-negative. .Sh RETURN VALUES These functions return the requested square root. @@ -83,12 +88,15 @@ an \*(Na is generated, an invalid exception will be thrown. .Xr math 3 , .Sh STANDARDS The -.Fn csqrt +.Fn csqrt , +.Fn csqrtf , and -.Fn csqrtf +.Fn csqrtl functions conform to .St -isoC-99 . .Sh BUGS For -.Fn csqrt , -inexact results are not correctly rounded in general. +.Fn csqrt +and +.Fn csqrtl , +inexact results are not always correctly rounded. diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c index 50e24ba..ae38c6c 100644 --- a/lib/msun/src/s_csqrt.c +++ b/lib/msun/src/s_csqrt.c @@ -28,6 +28,7 @@ __FBSDID("$FreeBSD$"); #include +#include #include #include "math_private.h" @@ -105,3 +106,7 @@ csqrt(double complex z) else return (result); } + +#if LDBL_MANT_DIG == 53 +__weak_reference(csqrt, csqrtl); +#endif diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c new file mode 100644 index 0000000..5a5087f --- /dev/null +++ b/lib/msun/src/s_csqrtl.c @@ -0,0 +1,108 @@ +/*- + * Copyright (c) 2007-2008 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 + +#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 >= LDBL_MAX / (1 + sqrt(2)). */ +#define THRESH (LDBL_MAX / 2.414213562373095048801688724209698L) + +long double complex +csqrtl(long double complex z) +{ + long double complex result; + long double a, b; + long double t; + int scale; + + a = creall(z); + b = cimagl(z); + + /* Handle special cases. */ + if (z == 0) + return (cpackl(0, b)); + if (isinf(b)) + return (cpackl(INFINITY, b)); + if (isnan(a)) { + t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ + return (cpackl(a, 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 (cpackl(fabsl(b - b), copysignl(a, b))); + else + return (cpackl(a, copysignl(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 (fabsl(a) >= THRESH || fabsl(b) >= THRESH) { + a *= 0.25; + b *= 0.25; + scale = 1; + } else { + scale = 0; + } + + /* Algorithm 312, CACM vol 10, Oct 1967. */ + if (a >= 0) { + t = sqrtl((a + hypotl(a, b)) * 0.5); + result = cpackl(t, b / (2 * t)); + } else { + t = sqrtl((-a + hypotl(a, b)) * 0.5); + result = cpackl(fabsl(b) / (2 * t), copysignl(t, b)); + } + + /* Rescale. */ + if (scale) + return (result * 2); + else + return (result); +} -- cgit v1.1