[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