Two possible problems with gsl_complex behavior

Two possible problems with gsl_complex behavior

[Charlie Zender]

Hi,

I use GSL on a RedHat GNU/Linux 7.0 system with GCC 2.96.
zender@cgd85:~/c++$ uname -a
Linux cgd85.cgd.ucar.edu 2.2.18 #2 Thu Mar 22 18:54:20 PST 2001 i686 unknown
zender@cgd85:~/c++$ g++ --version
2.96
I recently upgraded from gsl-0.7 to gsl-0.9.
There were two new behaviors that I noticed because they caused
a previously working program of mine to start failing. I'm not
saying these are bugs, perhaps they are features or compiler
problems, but, in any case, they appear to be worth investigating
before gsl-1.0. 

First, and least troublesome, is that the result of
gsl_complex_cos(gsl_complex z) where z=(0.0,0.0) does not appear to
be exactly (1.0,0.0), but rather (1.0,-0.0). 
In other words it prints as negative 0.0. 
Perhaps this is a feature, but to me it is a warning sign
that an argument to gsl_complex_cos that should give an exactly
representable answer does not.
The following test program illustrates this problem:

#include <iostream>
#include <complex>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

int main()
{
  gsl_complex ngl=gsl_complex_rect(0.0,0.0);
  gsl_complex ngl_cos=gsl_complex_cos(ngl);
  cout << "ngl = (" << GSL_REAL(ngl) << "," << GSL_IMAG(ngl) << ")" << endl;
  cout << "ngl_cos = (" << GSL_REAL(ngl_cos) << "," << GSL_IMAG(ngl_cos) << ")" << endl;
}

zender@cgd85:~/c++$ g++ -o tst tst.cc -L/usr/lib -lgsl -lgslcblas -lm
zender@cgd85:~/c++$ ./tst
ngl = (0,0)
ngl_cos = (1,-0)

Second, and more troublesome, is the behavior of gsl_complex functions
when used as arguments to other functions. Again, using g++-2.96, the
following used to work in gsl-0.7 "as expected":

float_complex cpx=(0.0,0.0);
gsl_complex ngl_cos=gsl_complex_cos(gsl_complex_rect(cpx.real(),cpx.imag())); 

Attempting this with gsl-0.9 causes an error, sometime a core dump and
sometimes a garbage result in ngl_cos. It does not appear to matter
whether I cast (using C++ style static_cast<double>() or C-style
(double)) the cpx.real() and cpx.imag() variables to doubles in the
arguments to gsl_complex_rect or not. For example, this causes a
core dump:

ngl_cos=gsl_complex_cos(gsl_complex_rect(static_cast<double>(cpx.real()),static_cast<double>(cpx.imag()))); 

Instead I must deconvolve the two functions as follows in order to
obtain "expected" behavior: 

float_complex cpx=(0.0,0.0);
gsl_complex ngl=gsl_complex_rect(cpx.real(),cpx.imag()); 
gsl_complex ngl_cos=gsl_complex_cos(ngl); 

In this deconvolved form it does not appear to be necessary to cast
cpx.real() and cpx.imag() to doubles before use in gsl_complex_rect().
Thus something appears to have changed regarding the intrinsic promotion
of arguments to the gsl_complex_rect() function. Of course
float_complex is based on floats, and gsl_complex is based on doubles
(did it used to be based on floats in gsl-0.7?), so I think what is
happening in gsl-0.9 is that the promotion of the float components of
cpx were not getting promoted to doubles as required by
gsl_complex_rect(). Any insight into this would be appreciated.

Overall, my experiences with GSL have been very positive.
Thanks for writing this software. I hope this report proves useful.

Charlie
-- 
Charlie Zender zender@uci.edu (949) 824-2987/FAX-3256, Department of
Earth System Science, University of California, Irvine CA 92697-3100
Visiting NCAR CGD Jul 1-Aug 31: Room ML-304A (303) 497-1738/FAX-1324
Work Address: NCAR CGD, 1850 Table Mesa Dr., Boulder, CO, 80303-5602

Re: Two possible problems with gsl_complex behavior

[Brian Gough]

Charlie Zender writes:
 > First, and least troublesome, is that the result of
 > gsl_complex_cos(gsl_complex z) where z=(0.0,0.0) does not appear to
 > be exactly (1.0,0.0), but rather (1.0,-0.0). 

Thanks for the bug reports.  

For cos(z) I've adjusted the code so that the sign of zero is positive
for Im(z) = 0.  

The negative sign was a result of the identity used to calculate the
imaginary part,

  cos(x+iy) = cosh(x) cos(y) - i sin(x) sinh(y)

 > Second, and more troublesome, is the behavior of gsl_complex functions
 > when used as arguments to other functions. Again, using g++-2.96, the
 > following used to work in gsl-0.7 "as expected":
 > 
 > float_complex cpx=(0.0,0.0);
 > gsl_complex ngl_cos=gsl_complex_cos(gsl_complex_rect(cpx.real(),cpx.imag())); 

I had a quick look and this function does not appear to have changed
between 0.7 and 0.9 apart from,

-- the name being changed from gsl_complex_xy() to gsl_complex_rect()

-- the addition of an inline version in the header file, wrapped with
#ifdef HAVE_INLINE (so this should not cause any change unless
compiling the test program with -DHAVE_INLINE)

I tried compiling with g++ 2.95 and it worked ok.

Before I investigate more did you try compiling with -Wall and other
warning options?  Also, you know that the initializer should be
float_complex cpx(0.0,0.0); rather than =(0.0,0.0); right?

regards
Brian Gough