Is cascading hypothesis () safe?

Question

I would like to calculate the norm (length) of three- and four-dimensional vectors. I use double precision floating point numbers and want to be careful to avoid unnecessary overflow or overflow.

The C math library provides hypot(x,y) for calculating the norm of two-dimensional vectors, trying to avoid overflow / overflow of intermediate calculations.

My question is: is it safe to use hypot(x, hypot(y, z)) and hypot(hypot(w, x), hypot(y, z)) to calculate the lengths of three- and four-dimensional vectors, respectively?

Answer 1

It is safe, but it is wasteful: you only need to compute sqrt() once, but when you cascade hypot() , you will call sqrt() for every call to hypot() . Usually, performance may not bother me, but it can also degrade the accuracy of the result. You can write your own:

 double hypot3(double x, double y, double z) { return sqrt(x*x + y*y + z*z); } 

etc .. It will be faster and more accurate. I don’t think anyone would be confused when they see hypot3() in your code.

The standard hypot() library may have tricks to avoid overflow, but you can not worry about that. Usually hypot() is more accurate than sqrt(x*x + y*y) . See e_hypot.c in the GLibC source code.

Answer 2

It is safe (almost) to use hypot(x, hypot(y, z)) and hypot(hypot(w, x), hypot(y, z)) to calculate the lengths of three- and four-dimensional vectors.

C does not indicate that hypot() should work for double x, y , which has the final answer of double . It has the words of affection "without excessive overflow or insufficient flow."

However, given that hypot(x, y) works, a reasonable implementation of hypot() will execute hypot(hypot(w, x), hypot(y, z)) as necessary. There is only 1 increment (at the lower end) / decrement (at the upper end) of the binary range of exponentials lost at 4-D and 2-D.

Regarding speed, accuracy and range, the code code from sqrtl((long double) w*w + (long double) x*x + (long double) y*y + (long double) z*z) is an alternative, but this It seems only necessary when choosing encoding targets.

Answer 3

I did some experiments with such things. In particular, I looked at a simple implementation, an implementation using hypots and (C-version of the reference version) of the BLAS DNRM2 function.

I found that in terms of overflow and downstream, the BLAS and hypot implementations were the same (in my tests) and far superior to the regular implementation. As for time, for large (hundreds) of vectors of size BLAS was about 6 times slower than the plain, while the hypothesis was 3 times slower than BLAS. Time differences were smaller for smaller sizes.

Answer 4

If the code cannot use hypot() or the wider precision types, the slow method checks the performance using frexp() and scales the @greggo arguments.

 #include <math.h> double nibot_norm(double w, double x, double y, double z) { // Sort the values by some means if (fabs(x) < fabs(w)) return nibot_norm(x, w, y, z); if (fabs(y) < fabs(x)) return nibot_norm(w, y, x, z); if (fabs(z) < fabs(y)) return nibot_norm(w, x, z, y); if (z == 0.0) return 0.0; // all zero case // Scale z to exponent half-way 1.0 to MAX_DOUBLE/4 // and w,x,y the same amount int maxi; frexp(DBL_MAX, &maxi); int zi; frexp(z, &zi); int pow2scale = (maxi / 2 - 2) - zi; // NO precision loss expected so far. // except w,x,y may become 0.0 if _far_ less than z w = ldexp(w, pow2scale); x = ldexp(x, pow2scale); y = ldexp(y, pow2scale); z = ldexp(z, pow2scale); // All finite values in range of squaring except for values // greatly insignificant to z (eg |z| > |x|*1e300) double norm = sqrt(((w * w + x * x) + y * y) + z * z); // Restore scale return ldexp(norm, -pow2scale); } 

Test code

 #include <float.h> #include <stdio.h> #ifndef DBL_TRUE_MIN #define DBL_TRUE_MIN DBL_MIN*DBL_EPSILON #endif void nibot_norm_test(double w, double x, double y, double z, double expect) { static int dig = DBL_DECIMAL_DIG - 1; printf(" w:%.*ex:%.*ey:%.*ez:%.*e\n", dig, w, dig, x, dig, y, dig, z); double norm = nibot_norm(w, x, y, z); printf("expect:%.*e\n", dig, expect); printf("actual:%.*e\n", dig, norm); if (expect != norm) puts("Different"); } int main(void) { nibot_norm_test(0, 0, 0, 0, 0); nibot_norm_test(10 / 7., 4 / 7., 2 / 7., 1 / 7., 11 / 7.); nibot_norm_test(DBL_MAX, 0, 0, 0, DBL_MAX); nibot_norm_test(DBL_MAX / 2, DBL_MAX / 2, DBL_MAX / 2, DBL_MAX / 2, DBL_MAX); nibot_norm_test(DBL_TRUE_MIN, 0, 0, 0, DBL_TRUE_MIN); nibot_norm_test(DBL_TRUE_MIN, DBL_TRUE_MIN, DBL_TRUE_MIN, DBL_TRUE_MIN, DBL_TRUE_MIN * 2); return 0; } 

results

  w:0.00000000000000000e+00 x:0.00000000000000000e+00 y:0.00000000000000000e+00 z:0.00000000000000000e+00 expect:0.00000000000000000e+00 actual:0.00000000000000000e+00 w:1.42857142857142860e+00 x:5.71428571428571397e-01 y:2.85714285714285698e-01 z:1.42857142857142849e-01 expect:1.57142857142857140e+00 actual:1.57142857142857140e+00 w:1.79769313486231571e+308 x:0.00000000000000000e+00 y:0.00000000000000000e+00 z:0.00000000000000000e+00 expect:1.79769313486231571e+308 actual:1.79769313486231571e+308 w:8.98846567431157854e+307 x:8.98846567431157854e+307 y:8.98846567431157854e+307 z:8.98846567431157854e+307 expect:1.79769313486231571e+308 actual:1.79769313486231571e+308 w:4.94065645841246544e-324 x:0.00000000000000000e+00 y:0.00000000000000000e+00 z:0.00000000000000000e+00 expect:4.94065645841246544e-324 actual:4.94065645841246544e-324 w:4.94065645841246544e-324 x:4.94065645841246544e-324 y:4.94065645841246544e-324 z:4.94065645841246544e-324 expect:9.88131291682493088e-324 actual:9.88131291682493088e-324