standards/104841: C99 long double square root.

Steven G. Kargl kargl at troutmask.apl.washington.edu
Thu Oct 26 21:30:21 UTC 2006


>Number:         104841
>Category:       standards
>Synopsis:       C99 long double square root.
>Confidential:   no
>Severity:       non-critical
>Priority:       medium
>Responsible:    freebsd-standards
>State:          open
>Quarter:        
>Keywords:       
>Date-Required:
>Class:          change-request
>Submitter-Id:   current-users
>Arrival-Date:   Thu Oct 26 21:30:17 GMT 2006
>Closed-Date:
>Last-Modified:
>Originator:     Steven G. Kargl
>Release:        FreeBSD 7.0-CURRENT amd64
>Organization:
APL/UW	
>Environment:
System: FreeBSD troutmask.apl.washington.edu 7.0-CURRENT FreeBSD 7.0-CURRENT #0: Tue Oct 17 13:17:15 PDT 2006 kargl at troutmask.apl.washington.edu:/usr/obj/usr/src/sys/SPEW amd64


	
>Description:

FreeBSD lacks an implementation of the C99 long double square root.

>How-To-Repeat:
	
>Fix:

The attach patch is an implementation of a long double square root.
The algorithm uses 2 iterations of Heron's algorithm after providing
an initial guess of the value from (almost) minmax polynomials.  These
polynomials provide an estimate with absolute error on the order of
1E-6.

An exhaustive evaluation of sqrtl() over all floats shows that the
ULP is less than or equal to 0.5.  An additional test of 10 million
long double values drawn from /dev/random also give <= 0.5 ULP.
These ULP tests compared sqrtl() to the value computed by GMP/MPFR
with 128 bits of precision, and then converted to a long double.


diff -ruN msun.orig/Makefile msun/Makefile
--- msun.orig/Makefile	Thu Apr 13 09:05:37 2006
+++ msun/Makefile	Thu Oct 26 14:03:27 2006
@@ -69,7 +69,7 @@
 .endif
 
 # C99 long double functions
-COMMON_SRCS+=	s_copysignl.c s_fabsl.c
+COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_sqrtl.c
 .if ${LDBL_PREC} != 53
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=	s_fmal.c s_frexpl.c s_nextafterl.c s_nexttoward.c s_scalbnl.c
diff -ruN msun.orig/man/sqrt.3 msun/man/sqrt.3
--- msun.orig/man/sqrt.3	Fri Jan 14 15:28:28 2005
+++ msun/man/sqrt.3	Thu Oct 26 14:01:18 2006
@@ -39,7 +39,8 @@
 .Nm cbrt ,
 .Nm cbrtf ,
 .Nm sqrt ,
-.Nm sqrtf
+.Nm sqrtf ,
+.Nm sqrtl
 .Nd cube root and square root functions
 .Sh LIBRARY
 .Lb libm
@@ -53,6 +54,8 @@
 .Fn sqrt "double x"
 .Ft float
 .Fn sqrtf "float x"
+.Ft long double
+.Fn sqrtl "long double x"
 .Sh DESCRIPTION
 The
 .Fn cbrt
@@ -63,9 +66,10 @@
 .Ar x .
 .Pp
 The
-.Fn sqrt
+.Fn sqrt ,
+.Fn sqrtf ,
 and the
-.Fn sqrtf
+.Fn sqrtl
 functions compute the
 non-negative square root of x.
 .Sh RETURN VALUES
@@ -75,9 +79,10 @@
 .Fn cbrtf
 functions return the requested cube root.
 The
-.Fn sqrt
+.Fn sqrt ,
+.Fn sqrtf ,
 and the
-.Fn sqrtf
+.Fn sqrtl
 functions return the requested square root
 unless an error occurs.
 An attempt to take the
diff -ruN msun.orig/src/math.h msun/src/math.h
--- msun.orig/src/math.h	Sat Apr 16 14:12:47 2005
+++ msun/src/math.h	Thu Oct 26 14:09:39 2006
@@ -460,7 +460,9 @@
 #if 0
 long double	sinhl(long double);
 long double	sinl(long double);
+#endif
 long double	sqrtl(long double);
+#if 0
 long double	tanhl(long double);
 long double	tanl(long double);
 long double	tgammal(long double);
diff -ruN msun.orig/src/s_sqrtl.c msun/src/s_sqrtl.c
--- msun.orig/src/s_sqrtl.c	Wed Dec 31 16:00:00 1969
+++ msun/src/s_sqrtl.c	Thu Oct 26 13:58:12 2006
@@ -0,0 +1,88 @@
+/*-
+ * Copyright (c) 2005, 2006 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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.
+ */
+
+/*
+ * Compute the sqrt root of a long double argument by decomposing the
+ * the argument into its base 2 fraction and exponent x = f 2**n where
+ * f = [1/2,1).
+ *
+ * If n is even, then sqrtl(x) = sqrt(f)   * sqrt(2**n)
+ * If n is odd,  then sqrtl(x) = sqrt(2*f) * sqrt(2**(n-1)) 
+ *
+ * The square root is computed via two iterations of Heron's algorithm
+ * where sqrt(f) and sqrt(2*f) are initially approximated by almost
+ * minmax polynomials with magnitudes of the absolute error on the
+ * order of 1E-6.
+ */
+
+#include <math.h>
+
+#define A0	 0x1.4cba1cbe49e33p-2
+#define A1	 0x1.d577769eba48bp-1
+#define A2	-0x1.46375f4632393p-2
+#define A3	 0x1.6551baf747915p-4
+#define A4	-0x1.5871c38c4c1ecp-7
+
+#define B0	 0x1.d68c0aa668916p-3
+#define B1	 0x1.4bf68abbcae36p+0
+#define B2	-0x1.cd56ea45d3571p-1
+#define B3	 0x1.f95363d9c4becp-2
+#define B4	-0x1.e71e308658219p-4
+
+long double
+sqrtl(long double x)
+{
+	int n;
+	long double f, s;
+
+	/* sqrt(NaN) = NaN, sqrt(+Inf) = +Inf, or  sqrt(-Inf) = sNaN. */
+	if (isnan(x) || isinf(x))
+		return (x * x + x);
+      
+	/* sqrt(+-0) = 0. */
+	if (x == 0.L)
+		return (0.L);
+
+	/* sqrt(x), x < 0. */
+	if (x < 0.L)
+		return ((x - x) / (x - x));
+
+	f = frexpl(x, &n);
+
+	if (n & 1) {
+		f *= 2;
+		n--;
+		s = A0 + (A1 + (A2 + (A3 + A4 * f) * f) * f) * f;
+	} else
+		s = B0 + (B1 + (B2 + (B3 + B4 * f) * f) * f) * f;
+
+	/* Heron's algorithm. */
+	s += f / s;
+	s = ldexpl(s, -1);
+	s += f / s;
+
+	return (ldexpl(s, (n >> 1) - 1));
+}
>Release-Note:
>Audit-Trail:
>Unformatted:


More information about the freebsd-standards mailing list