From owner-freebsd-standards@FreeBSD.ORG Fri Oct 12 18:14:27 2007 Return-Path: Delivered-To: freebsd-standards@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id ED43B16A417 for ; Fri, 12 Oct 2007 18:14:27 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.freebsd.org (Postfix) with ESMTP id C4FAC13C468 for ; Fri, 12 Oct 2007 18:14:27 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.1/8.14.1) with ESMTP id l9CI9xNN036377 for ; Fri, 12 Oct 2007 11:09:59 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.1/8.14.1/Submit) id l9CI9xR3036376 for freebsd-standards@freebsd.org; Fri, 12 Oct 2007 11:09:59 -0700 (PDT) (envelope-from sgk) Date: Fri, 12 Oct 2007 11:09:59 -0700 From: Steve Kargl To: freebsd-standards@freebsd.org Message-ID: <20071012180959.GA36345@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="VS++wcV0S1rZb1Fb" Content-Disposition: inline User-Agent: Mutt/1.4.2.3i Subject: [PATCH] hypotl, cabsl, and code removal in cabs X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 12 Oct 2007 18:14:28 -0000 --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, + * Placed into the Public Domain, 1994. + * + * Modified by Steven G. Kargl for the long double type. + */ + +#include +#include + +long double +cabsl(long double z) +{ + return hypotl(creall(z), cimagl(z)); +} --VS++wcV0S1rZb1Fb--