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