Why hypot() function is so slow?

I found Math.hypot() to be abysmally slow. I found I could code a fast java version using the same algorithm that produces identical results. This can be used to replace it.

/**
 * <b>hypot</b>
 * @param x
 * @param y
 * @return sqrt(x*x +y*y) without intermediate overflow or underflow. 
 * @Note {@link Math#hypot} is unnecessarily slow.  This returns the identical result to 
 * Math.hypot with reasonable run times (~40 nsec vs. 800 nsec). 
 * <p>The logic for computing z is copied from "Freely Distributable Math Library" 
 * fdlibm's e_hypot.c. This minimizes rounding error to provide 1 ulb accuracy.
 */
public static double hypot(double x, double y) {

    if (Double.isInfinite(x) || Double.isInfinite(y)) return Double.POSITIVE_INFINITY;
    if (Double.isNaN(x) || Double.isNaN(y)) return Double.NaN;

    x = Math.abs(x);
    y = Math.abs(y);

    if (x < y) {
        double d = x;
        x = y;
        y = d;
    }

    int xi = Math.getExponent(x);
    int yi = Math.getExponent(y);

    if (xi > yi + 27) return x;

    int bias = 0;
    if (xi > 510 || xi < -511) {
        bias = xi;
        x = Math.scalb(x, -bias);
        y = Math.scalb(y, -bias);           
    }


    // translated from "Freely Distributable Math Library" e_hypot.c to minimize rounding errors
    double z = 0; 
    if (x > 2*y) { 
        double x1 = Double.longBitsToDouble(Double.doubleToLongBits(x) & 0xffffffff00000000L);
        double x2 = x - x1;
        z = Math.sqrt(x1*x1 + (y*y + x2*(x+x1)));
    } else {
        double t = 2 * x;
        double t1 = Double.longBitsToDouble(Double.doubleToLongBits(t) & 0xffffffff00000000L);
        double t2 = t - t1;
        double y1 = Double.longBitsToDouble(Double.doubleToLongBits(y) & 0xffffffff00000000L);
        double y2 = y - y1;
        double x_y = x - y;
        z = Math.sqrt(t1*y1 + (x_y*x_y + (t1*y2 + t2*y))); // Note: 2*x*y <= x*x + y*y
    }

    if (bias == 0) {
        return z; 
    } else {
        return Math.scalb(z, bias);
    }
}

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 convergence

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

EDIT (Update from J.M):

Here is the original Moler-Morrison paper, and here is a nice follow-up due to Dubrulle.

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