[Bug-gsl] Replace gsl_hypot() with Moler & Morrison algortihm?

[James Ward]

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