Date: Fri, 12 Oct 2007 11:09:59 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-standards@freebsd.org Subject: [PATCH] hypotl, cabsl, and code removal in cabs Message-ID: <20071012180959.GA36345@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
--VS++wcV0S1rZb1Fb
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
The attached patch implements hypotl(), cabsl(), and removes
z_abs() from w_cabs.c. z_abs() is undocumented namespace
pollution in libm, and no header file under /usr/include
provides a prototype. The documentation has been updated,
and in particular an incorrect paragraph in cabs.3 has been
removed.
The patch does not contain the relevant diffs for Makefile,
Symbol.map, and math.h because it has become too difficult
to separate out the related parts from the unrelated parts
of other previously submitted patches.
AS with other patches, I've only tested on amd64, so these
should be include in libm only if LDBL_PREC == 64.
--
Steve
--VS++wcV0S1rZb1Fb
Content-Type: text/x-diff; charset=us-ascii
Content-Disposition: attachment; filename="msun1.diff"
Index: man/hypot.3
===================================================================
RCS file: /home/ncvs/src/lib/msun/man/hypot.3,v
retrieving revision 1.14
diff -u -p -r1.14 hypot.3
--- man/hypot.3 9 Jan 2007 01:02:06 -0000 1.14
+++ man/hypot.3 12 Oct 2007 17:58:39 -0000
@@ -34,8 +34,10 @@
.Sh NAME
.Nm hypot ,
.Nm hypotf ,
+.Nm hypotl ,
.Nm cabs ,
-.Nm cabsf
+.Nm cabsf ,
+.Nm cabsl
.Nd Euclidean distance and complex absolute value functions
.Sh LIBRARY
.Lb libm
@@ -45,25 +47,31 @@
.Fn hypot "double x" "double y"
.Ft float
.Fn hypotf "float x" "float y"
+.Ft "long double"
+.Fn hypotl "long double x" "long double y"
.In complex.h
.Ft double
.Fn cabs "double complex z"
.Ft float
.Fn cabsf "float complex z"
+.Ft "long double"
+.Fn cabsl "long double complex z"
.Sh DESCRIPTION
The
-.Fn hypot
-and
+.Fn hypot ,
.Fn hypotf
+and
+.Fn hypotl
functions
compute the
sqrt(x*x+y*y)
in such a way that underflow will not happen, and overflow
occurs only if the final result deserves it.
The
-.Fn cabs
-and
+.Fn cabs ,
.Fn cabsf
+and
+.Fn cabsl
functions compute the complex absolute value of
.Fa z .
.Pp
@@ -82,11 +90,6 @@ Consequently
exactly;
in general, hypot and cabs return an integer whenever an
integer might be expected.
-.Pp
-The same cannot be said for the shorter and faster version of hypot
-and cabs that is provided in the comments in cabs.c; its error can
-exceed 1.2
-.Em ulps .
.Sh NOTES
As might be expected,
.Fn hypot "v" "\*(Na"
Index: src/e_hypotl.c
===================================================================
RCS file: src/e_hypotl.c
diff -N src/e_hypotl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/e_hypotl.c 12 Oct 2007 17:58:39 -0000
@@ -0,0 +1,126 @@
+/*-
+ * 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.
+ */
+
+/*
+ * Compute hypotl(x, y) = sqrt(x**2 + y**2) such that underflow and
+ * overflow are avoided. To accomplishe the task, write x = a * 2**n
+ * and y = b * 2**m, one then has
+ *
+ * hypotl(x, y) = 2**n * sqrt(a**2 + b**(2*(m - n))) for n >= m
+ * or
+ * hypotl(x, y) = 2**m * sqrt(a**(2*(n - m)) + b**2) for n < m
+ *
+ * where a and b are in [0.5, 1). If (m - n) is less than -32 (for
+ * a 64-bit significand), then b**(2*(m - n)) can be ignored.
+ *
+ * Special cases:
+ * hypotl(+-Inf, NaN) = hypotl(NaN, +-Inf) = Inf
+ * hypotl(+-x, 0) = hypotl(0, +-y) = |x| + |y|
+ * hypotl(NaN, NaN) = hypotl(NaN, NaN) = Inf
+ */
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#define ZERO 0.L
+
+/*
+ * Changing this to LDBL_MANT_DIG / 2 allows hypotl to be used for other long
+ * double types with precision other than 64. This, of course, assumes a
+ * function sqrtl() function exists.
+ */
+#define PREC 32
+
+long double
+hypotl(long double x, long double y)
+{
+ union IEEEl2bits a, b;
+ int m, n;
+ long double val;
+
+ a.e = x;
+ a.bits.sign = 0;
+ b.e = y;
+ b.bits.sign = 0;
+
+ /* If x = 0 or y = 0, then hypotl(x, y) = |x| + |y|. */
+ if (!(a.bits.manh | a.bits.manl) || !(b.bits.manh | b.bits.manl))
+ return (a.e + b.e);
+ /*
+ * If x = +-Inf or y = +- Inf, then hypotl(x, y) = Inf (even if the
+ * other argument is NaN). If a.e = NaN or b.e = NaN and the other
+ * argument is not +- Inf, then hypotl(x, y) = NaN.
+ */
+ if (a.bits.exp == 32767 || b.bits.exp == 32767) {
+ mask_nbit_l(a);
+ mask_nbit_l(b);
+ if (!(a.bits.manh | a.bits.manl) ||
+ !(b.bits.manh | b.bits.manl))
+ return (1.L / ZERO);
+ return ((x - x) / (x - x));
+ }
+
+ /*
+ * At this point, a and b are normal or subnormal numbers and neither
+ * a nor b can be zero. Extract the exponents of a and b, and then
+ * reset the exponents such that a and b are in [0.5, 1).
+ */
+ if (a.bits.exp == 0) {
+ a.e *= 0x1.0p514;
+ n = a.bits.exp - 0x4200;
+ } else
+ n = a.bits.exp - 0x3ffe;
+ a.bits.exp = 0x3ffe;
+
+ if (b.bits.exp == 0) {
+ b.e *= 0x1.0p514;
+ m = b.bits.exp - 0x4200;
+ } else
+ m = b.bits.exp - 0x3ffe;
+ b.bits.exp = 0x3ffe;
+
+ if (n >= m) {
+ a.e *= a.e;
+ m -= n;
+ if (m > - PREC) {
+ b.bits.exp += m;
+ a.e += b.e * b.e;
+ }
+ a.e = sqrtl(a.e);
+ a.bits.exp += n;
+ return (a.e);
+ } else {
+ b.e *= b.e;
+ n -= m;
+ if (n > - PREC) {
+ a.bits.exp += n;
+ b.e += a.e * a.e;
+ }
+ b.e = sqrtl(b.e);
+ b.bits.exp += m;
+ return (b.e);
+ }
+}
Index: src/w_cabs.c
===================================================================
RCS file: /home/ncvs/src/lib/msun/src/w_cabs.c,v
retrieving revision 1.4
diff -u -p -r1.4 w_cabs.c
--- src/w_cabs.c 13 Jun 2001 15:16:30 -0000 1.4
+++ src/w_cabs.c 12 Oct 2007 17:58:39 -0000
@@ -19,10 +19,3 @@ cabs(z)
{
return hypot(creal(z), cimag(z));
}
-
-double
-z_abs(z)
- double complex *z;
-{
- return hypot(creal(*z), cimag(*z));
-}
Index: src/w_cabsl.c
===================================================================
RCS file: src/w_cabsl.c
diff -N src/w_cabsl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ src/w_cabsl.c 12 Oct 2007 17:58:39 -0000
@@ -0,0 +1,17 @@
+/*
+ * cabs() wrapper for hypot().
+ *
+ * Written by J.T. Conklin, <jtc@wimsey.com>
+ * Placed into the Public Domain, 1994.
+ *
+ * Modified by Steven G. Kargl for the long double type.
+ */
+
+#include <complex.h>
+#include <math.h>
+
+long double
+cabsl(long double z)
+{
+ return hypotl(creall(z), cimagl(z));
+}
--VS++wcV0S1rZb1Fb--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20071012180959.GA36345>