Date: Sat, 3 Nov 2007 16:15:27 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-standards@FreeBSD.org Subject: Re: [PATCH] hypotl, cabsl, and code removal in cabs Message-ID: <20071103155721.K671@besplex.bde.org> In-Reply-To: <20071012180959.GA36345@troutmask.apl.washington.edu> References: <20071012180959.GA36345@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20071103155721.K671>