[Bug-gsl] Complex: gsl_complex_div_real bug

[chris coleman-smith]

Hi
There is a bug in the implementation of gsl_complex_div_real  
(gsl_complex a, double x), from the 1.9 package in the complex/math.c  
file.

The code currently is:

gsl_complex
gsl_complex_div_real (gsl_complex a, double x)
{                               /* z=a/x */
  gsl_complex z;
  GSL_SET_COMPLEX (&z, GSL_REAL (a) / x, GSL_IMAG (a) / x);
  return z;
}


However in the situation where GSL_REAL(a) != 0 && GSL_IMAG(a) == 0  
and x != 0 we get the imag part of a set to INF when really it should  
be zero. Since the imag part of x is actually zero.

Perhaps this is simply a interpretation of the intent of this  
function, it does work if you're just trying to do a simple scalar  
divide of the complex number. I thought it was meant to provide the  
functionality of dividing a number in C by a number in R ( which is  
also in C) but it doesn't quite work.

Suggested change would be to either point out that this is intended to  
be a purely scalar divide (and not really reflective of how C / R  
division should work for which you should use the C/C division and  
cast your float into a C number) or put in some kind of simple if
        
GSL_IMAG(a) == 0 {
        GSL_SET_COMPLEX(&z, GSL_REAL(a)/x, 0);
}

Thanks

Chris