[PATCH] hypotl, cabsl, and code removal in cabs

[Steve Kargl]

On Sat, Nov 03, 2007 at 04:15:27PM +1100, Bruce Evans wrote:
> 
> 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.

Ah, no.  Underflow is prevented by the if() statement. 

		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);

If m < -PREC, then |b|/|a| < 2**(-32), so it skips the "a.e += b.e * b.e;"
"a.e. *= a.e" can't underflow because a.e is in [0.5,1), and 
"a.e = sqrtl(a.e);" can't underflow because, here, a.e is in [0.25,2).

> % 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 .

This paragraph is useless.  There is no cabs.c in the source tree,
and in particular, there is no cabs.c in an installed system.  Surely,
this paragraph isn't referring to src/contrib/libf2c/libF77/cabs.c.
It has one comment with a total of one word.

> 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.

I don't understand to what the above comment is directed.

> % 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).

See my comment above.  The usual case isn't the important case.
It's the unusual case that is of concern.  The above code prevents
the unusual, yet provides acceptable performance in the usual case.

> 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.

The above seems more complicated than needed when one considers
that x and y can be easily split into |x|=a*2**n and |y|=b*2**m
where f,g are in [0.5,1).   One then has the comment that is
found in the code:

 * Compute hypotl(x, y) = sqrt(x**2 + y**2) such that underflow and 
 * overflow are avoided.  To accomplish the task, write x = a * 2**n
 * and y = b * 2**m, one then has
 *
 *   hypotl(x, y) = 2**n * sqrt(a**2 + b**(2*(m - n))) for n >= m
 * or
 *   hypotl(x, y) = 2**m * sqrt(a**(2*(n - m)) + b**2) for n < m
 * 
 * where a and b are in [0.5, 1).  If (m - n) is less than -32 (for
 * a 64-bit significand), then b**(2*(m - n)) can be ignored.

-- 
Steve