views:

198

answers:

2

I did some testing with C++ hypot() and Java Math.hypot. They both seem to be significantly slower than sqrt(a*a + b*b). Is that because of a better precision? What method to calculate a hypotenuse hypot function uses? Surprisingly I couldn't find any indication of poor performance in the documentation.

+10  A: 

It's not a simple sqrt function. You should check this link for the implementation of the algorithm: http://www.koders.com/c/fid7D3C8841ADC384A5F8DE0D081C88331E3909BF3A.aspx

It has while loop to check for covergence

/* Slower but safer algorithm due to Moler and Morrison.  Never
         produces any intermediate result greater than roughly the
         larger of X and Y.  Should converge to machine-precision
         accuracy in 3 iterations.  */

      double r = ratio*ratio, t, s, p = abig, q = asmall;

      do {
        t = 4. + r;
        if (t == 4.)
          break;
        s = r / t;
        p += 2. * s * p;
        q *= s;
        r = (q / p) * (q / p);
      } while (1);
Ravi Gummadi
Would be nice to give examples where this function is better than the `sqrt` approach.
schnaader
@schnaader - read the linked page, the reason is right at the top (short version, the naive method may overflow when it shouldn't)
Nir
This matters for very large values of x and y. For example, if x and y are such that x^2 + y^2 is greater than MAX_DOUBLE, the sqrt(x^2 + y^2) version will fail. However, since this method never produces a value that is much larger than x or y, it should be safe in those situations.
PFHayes
Ah, thanks for clarification.
schnaader
A: 

Hello.

Here is a faster implementation, which results are also closer to java.lang.Math.hypot: (NB: for Delorie's implementation, need to add handling of NaN and +-Infinity inputs)

private static final double TWO_POW_450 = Double.longBitsToDouble(0x5C10000000000000L);
private static final double TWO_POW_N450 = Double.longBitsToDouble(0x23D0000000000000L);
private static final double TWO_POW_750 = Double.longBitsToDouble(0x6ED0000000000000L);
private static final double TWO_POW_N750 = Double.longBitsToDouble(0x1110000000000000L);
public static double hypot(double x, double y) {
    x = Math.abs(x);
    y = Math.abs(y);
    if (y < x) {
        double a = x;
        x = y;
        y = a;
    } else if (!(y >= x)) { // Testing if we have some NaN.
        if ((x == Double.POSITIVE_INFINITY) || (y == Double.POSITIVE_INFINITY)) {
            return Double.POSITIVE_INFINITY;
        } else {
            return Double.NaN;
        }
    }
    if (y-x == y) { // x too small to substract from y
        return y;
    } else {
        double factor;
        if (x > TWO_POW_450) { // 2^450 < x < y
            x *= TWO_POW_N750;
            y *= TWO_POW_N750;
            factor = TWO_POW_750;
        } else if (y < TWO_POW_N450) { // x < y < 2^-450
            x *= TWO_POW_750;
            y *= TWO_POW_750;
            factor = TWO_POW_N750;
        } else {
            factor = 1.0;
        }
        return factor * Math.sqrt(x*x+y*y);
    }
}
Jeff