Implement csqrtl().

5 files changed, 135 insertions, 13 deletions

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 <das@FreeBSD.org>
+.\" Copyright (c) 2007-2008 David Schultz <das@FreeBSD.org>
 .\" 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 <complex.h>
+#include <float.h>
 #include <math.h>
 
 #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 <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 <float.h>
+#include <math.h>
+
+#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);
+}