Date: Tue, 7 Aug 2007 17:58:55 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-current@freebsd.org Subject: Re: [PATCH] Annual sqrtl patch Message-ID: <20070808005855.GA9599@troutmask.apl.washington.edu> In-Reply-To: <20070806023232.GA47516@troutmask.apl.washington.edu> References: <20070806023232.GA47516@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
[-- Attachment #1 --]
On Sun, Aug 05, 2007 at 07:32:32PM -0700, Steve Kargl wrote:
> Here's the latest version of sqrtl. Please commit.
>
Here's a new patch that includes cbrtl, ccosh, and sqrtl.
--
Steve
[-- Attachment #2 --]
Index: include/complex.h
===================================================================
RCS file: /home/ncvs/src/include/complex.h,v
retrieving revision 1.6
diff -u -p -r1.6 complex.h
--- include/complex.h 14 Aug 2004 18:03:21 -0000 1.6
+++ include/complex.h 8 Aug 2007 00:56:41 -0000
@@ -45,6 +45,7 @@ __BEGIN_DECLS
double cabs(double complex);
float cabsf(float complex);
+double complex ccosh(double complex);
double cimag(double complex);
float cimagf(float complex);
long double cimagl(long double complex);
Index: lib/msun/Makefile
===================================================================
RCS file: /home/ncvs/src/lib/msun/Makefile,v
retrieving revision 1.78
diff -u -p -r1.78 Makefile
--- lib/msun/Makefile 21 May 2007 02:49:08 -0000 1.78
+++ lib/msun/Makefile 8 Aug 2007 00:56:59 -0000
@@ -68,13 +68,18 @@ 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 s_cbrtl.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
.endif
# C99 complex functions
-COMMON_SRCS+= s_cimag.c s_cimagf.c s_cimagl.c s_conj.c s_conjf.c s_conjl.c \
+COMMON_SRCS+= s_ccosh.c \
+ s_cimag.c s_cimagf.c s_cimagl.c s_conj.c s_conjf.c s_conjl.c \
s_creal.c s_crealf.c s_creall.c
# FreeBSD's C library supplies these functions:
Index: lib/msun/Symbol.map
===================================================================
RCS file: /home/ncvs/src/lib/msun/Symbol.map,v
retrieving revision 1.4
diff -u -p -r1.4 Symbol.map
--- lib/msun/Symbol.map 29 Apr 2007 14:05:20 -0000 1.4
+++ lib/msun/Symbol.map 8 Aug 2007 00:56:59 -0000
@@ -56,12 +56,15 @@ FBSD_1.0 {
sinhf;
sqrt;
sqrtf;
+ sqrtl;
asinh;
asinhf;
atan;
atanf;
cbrt;
cbrtf;
+ cbrtl;
+ ccosh;
ceil;
ceilf;
ceill;
Index: lib/msun/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
--- lib/msun/man/sqrt.3 9 Jan 2007 01:02:06 -0000 1.13
+++ lib/msun/man/sqrt.3 8 Aug 2007 00:56:59 -0000
@@ -28,14 +28,16 @@
.\" 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 cbrtl ,
.Nm sqrt ,
-.Nm sqrtf
+.Nm sqrtf ,
+.Nm sqrtl
.Nd cube root and square root functions
.Sh LIBRARY
.Lb libm
@@ -45,35 +47,43 @@
.Fn cbrt "double x"
.Ft float
.Fn cbrtf "float x"
+.Ft "long double"
+.Fn cbrtl "long double x"
.Ft double
.Fn sqrt "double x"
.Ft float
.Fn sqrtf "float x"
+.Ft "long double"
+.Fn sqrtl "long double x"
.Sh DESCRIPTION
The
-.Fn cbrt
-and the
-.Fn cbrtf
-functions compute
-the cube root of
+.Fn cbrt ,
+.Fn cbrtf ,
+and
+.Fn cbrtl
+functions compute the cube root of
.Ar x .
.Pp
The
.Fn sqrt
-and the
.Fn sqrtf
-functions compute the
-non-negative square root of x.
+and
+.Fn sqrtl
+functions compute the non-negative square root of x except
+.Fn sqrt "-0"
+= -0.
.Sh RETURN VALUES
The
-.Fn cbrt
+.Fn cbrt ,
+.Fn cbrtf ,
and the
-.Fn cbrtf
+.Fn cbrtl
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
Index: lib/msun/src/e_sqrtl.c
===================================================================
RCS file: lib/msun/src/e_sqrtl.c
diff -N lib/msun/src/e_sqrtl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/msun/src/e_sqrtl.c 8 Aug 2007 00:56:59 -0000
@@ -0,0 +1,92 @@
+#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;
+ fexcept_t flagp;
+
+ r = fegetround();
+ fegetexceptflag(&flagp, FE_INEXACT);
+
+ u.e = 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 x = +-0, then sqrt(x) = +-0. */
+ if ((u.bits.manh | u.bits.manl) == 0)
+ return (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. */
+ y = x / u.e; /* Chopped quotient (inexact?). */
+
+ if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */
+ if (y == u.e) {
+ /* Restore inexact and rounding mode. */
+ fesetexceptflag(&flagp, FE_INEXACT);
+ fesetround(r);
+ return (u.e);
+ }
+ y = nexttowardl(y, 0.); /* y = y - ulp. */
+ }
+
+ fegetexceptflag(&flagp, FE_INEXACT); /* sqrt(x) is inexact. */
+ if (r == FE_TONEAREST)
+ y = nexttowardl(y, x); /* y = y + ulp. */
+ if (r == FE_UPWARD) {
+ u.e = nexttowardl(u.e, x); /* u.e = u.e + ulp. */
+ y = nexttowardl(y, x); /* y = y + ulp. */
+ }
+ u.e = u.e + y; /* Chopped sum. */
+ u.bits.exp--; /* Correctly rounded. */
+ /* Restore inexact and rounding mode. */
+ fesetexceptflag(&flagp, FE_INEXACT);
+ fesetround(r);
+ return (u.e);
+}
Index: lib/msun/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
--- lib/msun/src/math.h 7 Jan 2007 07:54:21 -0000 1.62
+++ lib/msun/src/math.h 8 Aug 2007 00:56:59 -0000
@@ -397,8 +397,8 @@ long double asinl(long double);
long double atan2l(long double, long double);
long double atanhl(long double);
long double atanl(long double);
-long double cbrtl(long double);
#endif
+long double cbrtl(long double);
long double ceill(long double);
long double copysignl(long double, long double) __pure2;
#if 0
@@ -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);
Index: lib/msun/src/s_cbrtl.c
===================================================================
RCS file: lib/msun/src/s_cbrtl.c
diff -N lib/msun/src/s_cbrtl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/msun/src/s_cbrtl.c 8 Aug 2007 00:56:59 -0000
@@ -0,0 +1,90 @@
+/*-
+ * 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 "fpmath.h"
+
+long double
+cbrtl(long double x)
+{
+ union IEEEl2bits u;
+ int k;
+ long double e = 1., y, z;
+
+ u.e = x;
+
+ /* If x = +-Inf, then cbrt(x) = +-Inf. */
+ /* If x = NaN, then sqrt(x) = NaN. */
+ if (u.bits.exp == 32767)
+ return (x + x);
+
+ /* If x = +-0, then sqrt(x) = +-0. */
+ if ((u.bits.manh | u.bits.manl) == 0)
+ return (x);
+
+ if (u.bits.exp == 0) {
+ /* Adjust subnormal numbers. */
+ u.e *= 0x1.0p514;
+ k = u.bits.exp - 0x4200;
+ } else /* x is a normal number, remove the bias. */
+ k = u.bits.exp - 0x3ffe;
+
+ u.bits.exp = 0x3ffe;
+
+ /* x = u.e * 2^n. */
+ switch (k - 3 * (k / 3)) {
+ case 0: /* mod(n,3) == 0. */
+ k -= 3; /* k = n - 3. */
+ u.bits.exp += 3; /* u.e in [4,8). */
+ break;
+ case -2:
+ case 1: /* mod(n,3) == 1. */
+ k--; /* k = n - 1. */
+ u.bits.exp++; /* u.e in [1,2). */
+ break;
+ case -1:
+ case 2: /* mod(n,3) == 2. */
+ k -= 2; /* k = n - 2. */
+ u.bits.exp +=2; /* u.e in [2,4). */
+ break;
+ }
+ /*
+ * Newton's iteration. Split u.e into a high and low part
+ * to achieve additional precision.
+ */
+ y = cbrt(u.e);
+ z = 1 / (y * y);
+ e = u.e;
+ u.bits.manl = 0;
+ e = y + (e - u.e) * z;
+ y += u.e * z;
+ u.e = (y + e) / 3;
+ u.bits.exp += (k / 3);
+
+ return (u.e);
+}
Index: lib/msun/src/s_ccosh.c
===================================================================
RCS file: lib/msun/src/s_ccosh.c
diff -N lib/msun/src/s_ccosh.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/msun/src/s_ccosh.c 8 Aug 2007 00:56:59 -0000
@@ -0,0 +1,127 @@
+/*-
+ * Copyright (c) 2005 Bruce D. Evans and 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.
+ */
+
+/*
+ * Hyperbolic cosine of a complex argument z = x + i y.
+ *
+ * cosh(z) = cosh(x+iy)
+ * = cosh(x) cos(y) + i sinh(x) sin(y).
+ *
+ * Exceptional values are noted in the comments within the source code.
+ * These values and the return value were taken from n1124.pdf.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+#include "math_private.h"
+
+double complex
+ccosh(double complex z)
+{
+ double x, y;
+ int32_t hx, hy, ix, iy, lx, ly;
+
+ x = creal(z);
+ y = cimag(z);
+
+ EXTRACT_WORDS(hx, lx, x);
+ EXTRACT_WORDS(hy, ly, y);
+
+ ix = 0x7fffffff & hx; /* If ix >= 0x7ff00000, then inf or NaN. */
+ iy = 0x7fffffff & hy; /* If iy >= 0x7ff00000, then inf or NaN. */
+
+ /* Handle the nearly-non-exceptional cases where x and y are finite. */
+ if (ix < 0x7ff00000 && iy < 0x7ff00000) {
+ if ((iy | ly) == 0)
+ return (cpack(cosh(x), x * y));
+ return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
+ }
+
+ /*
+ * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
+ * The sign of 0 in the result is unspecified. Choice = normally
+ * the same as dNaN. Raise the invalid floating-point exception.
+ *
+ * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
+ * The sign of 0 in the result is unspecified. Choice = normally
+ * the same as d(NaN).
+ *
+ * Say something about our sign choices generally but not specifically?
+ * The above is specific; dNaN = default NaN and d(NaN) = input NaN
+ * with default transformation under operations (e.g., on i386's,
+ * convert signaling NaNs to quiet NaNs by setting the quiet bit).
+ */
+ if ((ix | lx) == 0 && iy >= 0x7ff00000)
+ return (cpack(y - y, copysign(0, x * (y - y))));
+
+ /*
+ * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
+ *
+ * cosh(NaN + I 0) = NaN +- I 0.
+ * The sign of 0 in the result is unspecified.
+ */
+ if ((iy | ly) == 0 && ix >= 0x7ff00000) {
+ if (((hx & 0xfffff) | lx) == 0)
+ return (cpack(x * x, copysign(0, x) * y));
+ return (cpack(x * x, copysign(0, (x + x) * y)));
+ }
+
+ /*
+ * cosh(x +- I Inf) = dNaN + I dNaN.
+ * Raise the invalid floating-point exception for finite nonzero x.
+ *
+ * cosh(x + I NaN) = d(NaN) + I d(NaN).
+ * Optionally raises the invalid floating-point exception for finite
+ * nonzero x. Choice = raise.
+ */
+ if (ix < 0x7ff00000 && iy >= 0x7ff00000)
+ return (cpack(y - y, x * (y - y)));
+
+ /*
+ * cosh(+-Inf + I NaN) = +Inf + I d(NaN).
+ *
+ * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
+ * The sign of Inf in the result is unspecified. Choice = always +.
+ * Raise the invalid floating-point exception.
+ *
+ * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
+ */
+ if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
+ if (iy >= 0x7ff00000)
+ return (cpack(x * x, x * (y - y)));
+ return (cpack((x * x) * cos(y), x * sin(y)));
+ }
+
+ /*
+ * cosh(NaN + I NaN) = NaN + I NaN.
+ *
+ * cosh(NaN + I y) = NaN + I NaN.
+ * Raise the invalid floating-point exception for finite nonzero y.
+ *
+ */
+ return (cpack((x * x) * (y - y), (x + x) * (y - y)));
+}
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20070808005855.GA9599>