[Rd] hypot(x,y) instead of pythag(a,b) ?!

[Martin Maechler]

Some of you may have seen the  pythag() part in the R API definition in 
"Writing R Extensions" (source = doc/manual/R-exts.texinfo).

or followed the report and Prof. Brian Ripley's answer about pythag()'s
availability from R's binary.

As we say in above manual

  >>      `pythag(A, B)' computes `sqrt(A^2 + B^2)' without overflow or
  >>      destructive underflow: for example it still works when both A and
  >>      B are between `1e200' and `1e300' (in IEEE double precision).

--
"Problem" is : 
The GNU C library (and other C libraries ??)
defines a function	 
	 double hypot(double x, double y)

with identical semantics to our pythag() from above
The Info (e.g. in Linux Emacs C-h i "m libc") about "Libc" contains
(in the section "Exponentiation and Logarithms"):

=============================
>>  - Function: double hypot (double X, double Y)
>>  - Function: float hypotf (float X, float Y)
>>  - Function: long double hypotl (long double X, long double Y)
>>      These functions return `sqrt (X*X + Y*Y)'.  This is the length of
>>      the hypotenuse of a right triangle with sides of length X and Y,
>>      or the distance of the point (X, Y) from the origin.  Using this
>>      function instead of the direct formula is wise, since the error is
>>      much smaller.  See also the function `cabs' in *Note Absolute
>>      Value::.

Further "problem": In R, we are already partially relying on the
availability of the hypot() function :

At the toplevel of R-1.0.1's source
   grep -rwn hypot . 
   ~~~~~~~~~~~~~~~~~ (with a newer GNU grep that has "-r" for "recursive"):
gives

./src/appl/cpoly.c:145:	shr[i] = hypot(pr[i], pi[i]);
./src/appl/fortran.c:111:    return hypot(z->r, z->i);
./src/main/complex.c:122:    logr = log(hypot(a->r, a->i) );
./src/main/complex.c:279:  REAL(y)[i] = hypot(COMPLEX(x)[i].r, COMPLEX(x)[i].i);
./src/main/complex.c:285:  REAL(y)[i] = hypot(COMPLEX(x)[i].r, COMPLEX(x)[i].i);
./src/main/complex.c:388:    r->r = log(hypot( z->r, z->i ));
./src/main/complex.c:411:    if( (mag = hypot(z->r, z->i)) == 0.0)
./src/main/plot.c:1201:		if ((f = d/hypot(xx-xold, yy-yold)) < 0.5) {
./src/main/plot.c:2455:double hypot(double x, double y)
./src/main/plot.c:2559:		d = hypot(xp-xi, yp-yi);
./src/gnuwin32/math/protos.h:43:extern double hypot ( double x, double y );

---------

My "theses"
   
   o when hypot() is available, it should be used since it will be 
     efficient, precise, etc.
   o when it's not available, our "configure" should find out, and set
     corresponding "HAVE_HYPOT" to `false'.
   o in that case, we should provide a hypot()  {for the above code to
     work!}, and we probably should use (an improvement?) of the current
     pythag() function [src/appl/pythag.c in R versions <= 1.0.1].
   
   o Hence, we should drop pythag() from the R API and rather say that we
     provide hypot() whenever the system's C library 
     (or its "math.h", aka libm.*, aka "-m" part, respectively) does *not* 
     provide it.

   o I think an R function  hypot() would also make sense.

We'll be glad to hear feedback about this!

Martin Maechler <maechler@stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO D10	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._