[PATCH] hypotl, cabsl, and code removal in cabs
Bruce Evans
brde at optusnet.com.au
Fri Nov 2 22:15:46 PDT 2007
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
More information about the freebsd-standards
mailing list