Date: Sun, 5 Aug 2007 19:32:32 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-current@freebsd.org Subject: [PATCH] Annual sqrtl patch Message-ID: <20070806023232.GA47516@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
--PNTmBPCT7hxwcZjr
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
Here's the latest version of sqrtl. Please commit.
--
Steve
--PNTmBPCT7hxwcZjr
Content-Type: text/x-diff; charset=us-ascii
Content-Disposition: attachment; filename="sqrtl.diff"
? sqrtl.diff
Index: Makefile
===================================================================
RCS file: /home/ncvs/src/lib/msun/Makefile,v
retrieving revision 1.78
diff -u -p -r1.78 Makefile
--- Makefile 21 May 2007 02:49:08 -0000 1.78
+++ Makefile 6 Aug 2007 02:30:58 -0000
@@ -68,6 +68,10 @@ SYMBOL_MAPS= ${SYM_MAPS}
# C99 long double functions
COMMON_SRCS+= s_copysignl.c s_fabsl.c s_modfl.c
+.if ${MACHINE_ARCH} == "i386" || ${MACHINE_ARCH} == "amd64"
+COMMON_SRCS+= e_sqrtl.c
+.endif
+
.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
Index: Symbol.map
===================================================================
RCS file: /home/ncvs/src/lib/msun/Symbol.map,v
retrieving revision 1.4
diff -u -p -r1.4 Symbol.map
--- Symbol.map 29 Apr 2007 14:05:20 -0000 1.4
+++ Symbol.map 6 Aug 2007 02:30:58 -0000
@@ -56,6 +56,7 @@ FBSD_1.0 {
sinhf;
sqrt;
sqrtf;
+ sqrtl;
asinh;
asinhf;
atan;
Index: man/sqrt.3
===================================================================
RCS file: /home/ncvs/src/lib/msun/man/sqrt.3,v
retrieving revision 1.13
diff -u -p -r1.13 sqrt.3
--- man/sqrt.3 9 Jan 2007 01:02:06 -0000 1.13
+++ man/sqrt.3 6 Aug 2007 02:30:58 -0000
@@ -28,14 +28,15 @@
.\" from: @(#)sqrt.3 6.4 (Berkeley) 5/6/91
.\" $FreeBSD: src/lib/msun/man/sqrt.3,v 1.13 2007/01/09 01:02:06 imp Exp $
.\"
-.Dd May 6, 1991
+.Dd August 5, 2007
.Dt SQRT 3
.Os
.Sh NAME
.Nm cbrt ,
.Nm cbrtf ,
.Nm sqrt ,
-.Nm sqrtf
+.Nm sqrtf ,
+.Nm sqrtl
.Nd cube root and square root functions
.Sh LIBRARY
.Lb libm
@@ -50,6 +51,9 @@
.Ft float
.Fn sqrtf "float x"
.Sh DESCRIPTION
+.Ft "long double"
+.Fn sqrtl "long double x"
+.Sh DESCRIPTION
The
.Fn cbrt
and the
Index: src/e_sqrtl.c
===================================================================
RCS file: src/e_sqrtl.c
diff -N src/e_sqrtl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/e_sqrtl.c 6 Aug 2007 02:30:58 -0000
@@ -0,0 +1,118 @@
+/*-
+ * Copyright (c) 2007, 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.
+ */
+
+#include <math.h>
+#include <stdint.h>
+
+#include <fenv.h>
+#pragma STDC FENV_ACCESS ON
+
+#include "fpmath.h"
+
+long double
+sqrtl(long double x)
+{
+ union IEEEl2bits u;
+ int k, r, i;
+ long double e, y, z;
+ fexcept_t flagp;
+
+ r = fegetround();
+ fegetexceptflag(&flagp, FE_INEXACT);
+
+ u.e = x;
+
+ /* If x = +-0, then sqrt(x) = +-0. */
+ if ((u.bits.manh | u.bits.manl) == 0)
+ return (x);
+
+ /* If x < 0, then sqrt(x) = sNaN */
+ if (u.bits.sign)
+ return ((x - x) / (x - x));
+
+ /* If x = NaN, then sqrt(x) = NaN. */
+ /* If x = Inf, then sqrt(x) = Inf. */
+ if (u.bits.exp == 32767)
+ return (x * x + x);
+
+ if (u.bits.exp == 0) {
+ /* Adjust subnormal numbers. */
+ u.e *= 0x1.0p514;
+ k = u.bits.exp - 0x4200;
+ u.bits.exp = 0x3ffe;
+ } else {
+ /*
+ * x is a normal number, so break it into x = e*2^n where
+ * x = (2*e)*2^2k for odd n and x = (4*e)*2^2k for even n.
+ */
+ if ((u.bits.exp - 0x3ffe) & 1) {/* n is odd. */
+ k = u.bits.exp - 0x3fff; /* 2k = n - 1. */
+ u.bits.exp = 0x3fff; /* u.e in [1,2). */
+ } else {
+ k = u.bits.exp - 0x4000; /* 2k = n - 2. */
+ u.bits.exp = 0x4000; /* u.e in [2,4). */
+ }
+ }
+ /*
+ * Newton's iteration. Split u.e into a high and low part
+ * to achieve additional precision.
+ */
+ y = sqrt(u.e);
+ e = u.e;
+ u.bits.manl = 0;
+ e = (e - u.e) / y;
+ y = (y + u.e / y);
+ u.e = y + e;
+ u.bits.exp += (k >> 1) - 1;
+
+ feclearexcept(FE_INEXACT); /* Reset INEXACT flag. */
+ fesetround(FE_TOWARDZERO); /* Set to round-toward-zero. */
+ z = x / u.e; /* Chopped quotient (inexact?). */
+
+ if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */
+ if (z == u.e) {
+ /* Restore inexact and rounding mode. */
+ fesetexceptflag(&flagp, FE_INEXACT);
+ fesetround(r);
+ return (u.e);
+ }
+ z = nexttowardl(z, 0.); /* z = z - ulp. */
+ }
+
+ fegetexceptflag(&flagp, FE_INEXACT); /* sqrt(x) is inexact. */
+ if (r == FE_TONEAREST)
+ z = nexttowardl(z, x); /* z = z + ulp. */
+ if (r == FE_UPWARD) {
+ u.e = nexttowardl(u.e, x); /* u.e = u.e + ulp. */
+ z = nexttowardl(z, x); /* z = z + ulp. */
+ }
+ u.e = u.e + z; /* Chopped sum. */
+ u.bits.exp--; /* Correctly rounded. */
+ /* Restore inexact and rounding mode. */
+ fesetexceptflag(&flagp, FE_INEXACT);
+ fesetround(r);
+ return (u.e);
+}
Index: src/math.h
===================================================================
RCS file: /home/ncvs/src/lib/msun/src/math.h,v
retrieving revision 1.62
diff -u -p -r1.62 math.h
--- src/math.h 7 Jan 2007 07:54:21 -0000 1.62
+++ src/math.h 6 Aug 2007 02:30:58 -0000
@@ -460,7 +460,9 @@ long double scalbnl(long double, int);
#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);
--PNTmBPCT7hxwcZjr--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20070806023232.GA47516>