[Rd] hypot(x,y) instead of pythag(a,b) ?!
Martin Maechler
Martin Maechler <maechler@stat.math.ethz.ch>
Wed, 24 May 2000 14:30:24 +0200 (CEST)
>>>>> "KH" == Kurt Hornik <Kurt.Hornik@ci.tuwien.ac.at> writes:
>>>>> Martin Maechler writes:
>> 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!
KH> I have no strong opinions about this, but the above seems to make
KH> sense.
KH> If I should make changes to configure, then pls let me know.
Yes, please.
We should check for it's presence with configure and
then use
#ifndef HAVE_HYPOT
double hypot(double x, double y) {
...
...
}
#endif
In the mean time I've seen that it exists in
- GNU lib"c" (libm infact)
- Solaris
- HP-UX
- DEC Alpha Unix
and it is part of the X/Open "XPG4" standard,
but not POSIX (as of Don Lewine) nor ANSI C.
Martin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._