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>
index | next in thread | previous in thread | raw e-mail
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
home |
help
Want to link to this message? Use this
URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20071103155721.K671>