Почему функция hypot() работает так медленно?

Question

Я провел некоторые испытания с помощью C++hypot() и JavaMath.hypot. Они кажутся значительно медленнее, чем sqrt(a*a + b*b). Это из-за лучшей точности? Какой метод расчета функции гипотенузы hypot использует? Удивительно, но я не мог найти никаких признаков плохой работы в документации.

Answer 1

Это не простая функция sqrt. Вы должны проверить эту ссылку для реализации алгоритма: http://www.koders.com/c/fid7D3C8841ADC384A5F8DE0D081C88331E3909BF3A.aspx

имеет время цикла для проверки сходимости

/* 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 (Обновление от JM):

Here является оригинальным Moler-Morrison бумага и here - прекрасное продолжение благодаря Dubrulle.

Answer 2

Вот быстрее реализация, результаты которого также ближе к java.lang.Math.hypot: (NB: для реализации Delorie, в нужно добавить обработку Нан и + -Infinity входы)

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

Answer 3

http://steve.hollasch.net/cgindex/math/pythag-root.txt

предполагает, что чем быстрее реализация с использованием SQRT() является квадратичной по сравнению с кубической для Moler & Morrison, примерно с теми же характеристиками переполнения.

Answer 4

Я нашел Math.hypot(), чтобы быть ужасно медленным. Я обнаружил, что могу скопировать быструю версию Java, используя тот же алгоритм, который дает идентичные результаты. Это можно использовать для его замены.

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