From owner-freebsd-standards@FreeBSD.ORG Sat Nov 3 05:15:46 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 4992A16A418 for ; Sat, 3 Nov 2007 05:15:46 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail35.syd.optusnet.com.au (mail35.syd.optusnet.com.au [211.29.133.51]) by mx1.freebsd.org (Postfix) with ESMTP id E1DF413C491 for ; Sat, 3 Nov 2007 05:15:45 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c211-30-219-213.carlnfd3.nsw.optusnet.com.au [211.30.219.213]) by mail35.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id lA35FRMG008246 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Sat, 3 Nov 2007 16:15:28 +1100 Date: Sat, 3 Nov 2007 16:15:27 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl In-Reply-To: <20071012180959.GA36345@troutmask.apl.washington.edu> Message-ID: <20071103155721.K671@besplex.bde.org> References: <20071012180959.GA36345@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@FreeBSD.org Subject: Re: [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: Sat, 03 Nov 2007 05:15:46 -0000 On Fri, 12 Oct 2007, Steve Kargl wrote: > 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. I fear that that paragraph is correct again, since the long double version only seems to take enough care with the overflow cases. % Index: man/hypot.3 % ... % @@ -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 . It's silly that on i386's, naive double precision code would work if the default precision were extended, but the double precision code isn't naive; OTOH, long double precision code can never depend on extra precision, but is naive. % Index: src/e_hypotl.c % ... % + 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); The usual case reduces to just sqrtl(x*x + y*y). fdlibm is more careful in this case (but similar in others): % * So, compute sqrt(x*x+y*y) with some care as % * follows to get the error below 1 ulp: % * % * Assume x>y>0; % * (if possible, set rounding to round-to-nearest) % * 1. if x > 2y use % * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y % * where x1 = x with lower 32 bits cleared, x2 = x-x1; else % * 2. if x <= 2y use % * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) % * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, % * y1= y with lower 32 bits chopped, y2 = y-y1. I don't understand the details of this, and couldn't find any cases where this makes a difference. Testing of 2-arg functions is hard since exhaustive testing is impossible even for floats, and the cases where this matters for hypot() seem to be sparse. Bruce