You could replace gsl_hypot() and gsl_hypot3() with Moler & Morrison based algorithms:
http://portal.acm.org/citation.cfm?id=1665041&CFID=15885790&CFTOKEN=32422364. But it might be slower than the current implementation. If you do decide to use M&M, you might want to fix anything that computes a Euclidean distance, like complex modulus or vector norm, to call gsl_hypot() instead of hypot(), which probably does what gsl_hypot() does now. I noticed that the current gsl_hypot() is not BSD on Infinity, NaN comparisons. The BSD manpage says that "hypot(+-infinity, NaN) = +infinity." http://www.ipnom.com/FreeBSD-Man-Pages/hypot.3.html but gsl_hypot() does not do this. If you add these 4 tests to ~/sys/test.c /* Test +Inf, NaN */ y = gsl_hypot (GSL_POSINF, GSL_NAN); y_expected = GSL_POSINF; gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_POSINF, GSL_NAN)"); /* Test -Inf, NaN */ y = gsl_hypot (GSL_NEGINF, GSL_NAN); y_expected = GSL_POSINF; gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NEGINF, GSL_NAN)"); /* Test NaN, +Inf */ y = gsl_hypot (GSL_NAN, GSL_POSINF); y_expected = GSL_POSINF; gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NAN, GSL_POSINF)"); /* Test NaN, -Inf */ y = gsl_hypot (GSL_NAN, GSL_NEGINF); y_expected = GSL_POSINF; gsl_test_rel (y, y_expected, 1e-15, "gsl_hypot(GSL_NAN, GSL_NEGINF)"); you will get: FAIL: gsl_hypot(GSL_POSINF, GSL_NAN) (nan observed vs inf expected) [27] FAIL: gsl_hypot(GSL_NEGINF, GSL_NAN) (nan observed vs inf expected) [29] FAIL: gsl_hypot(GSL_NAN, GSL_POSINF) (nan observed vs inf expected) [31] FAIL: gsl_hypot(GSL_NAN, GSL_NEGINF) (nan observed vs inf expected) [33] I found a C standard which says this Infinity, NaN behavior is optional: http://pubs.opengroup.org/onlinepubs/009695399/functions/hypot.html I guess you could: 1. Add a line to the Reference Manual saying that gsl_hypot() is not BSD on Infinity, NaN comparisons. or 2. Fix gsl_hypot() to be BSD. If you could get the BSD source code, you could see how they do it and add it to the current gsl_hypot(). Best wishes, Jim Ward _______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
