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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> "MM" == Martin Maechler <maechler@stat.math.ethz.ch> writes:MM> Some of you may have seen the pythag() part in the R API definition MM> in "Writing R Extensions" (source = doc/manual/R-exts.texinfo). MM> or followed the report and Prof. Brian Ripley's answer about MM> pythag()'s availability from R's binary. MM> 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). MM> -- "Problem" is : The GNU C library (and other C libraries ??) MM> defines a function double hypot(double x, double y) MM> with identical semantics to our pythag() from above The Info MM> (e.g. in Linux Emacs C-h i "m libc") about "Libc" contains (in the MM> section "Exponentiation and Logarithms"): .... I forgot to mention that src/appl/eigen.f (stemming from EISPACK) defines (and uses) a fortran DOUBLE PRECISION PYTHAG(A, B) function equivalently to src/appl/pythag.c {currently (R-devel) in .../approx.c}, and hence R's eigen() function actually does rely on pythag(). 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> 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!I have no strong opinions about this, but the above seems to make sense. If I should make changes to configure, then pls let me know. -k -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._