Re: 2.9.11?

[Daniel J Sebald]

Daniel J Sebald wrote:
> John W. Eaton wrote:
>
> > On 19-Apr-2007, Daniel J Sebald wrote:
> >
> > | So, it looks like the failing case is sqrt(-Inf).  OK, now where is 
> > the complex | version of square root, std::sqrt?  :-)
> >
> > It depends on your system.  With GCC, it would be in the complex
> > header file that is part of libstdc++.  Depending on compiler version
> > and options it could actually be defined using one of the C99
> > functions which are apparently built-in functions in GCC.
> >
> > Anyway, I don't think it is a bug in Octave.
>
>
> Oh, I didn't know there were complex versions of routines in the library:
>
> tgmath.h:#define sqrt(Val) __TGMATH_UNARY_REAL_IMAG (Val, sqrt, csqrt)
>
> So, it is really the csqrt that is being called, I guess.  Not the 
> sqrt() which would give an error for negative numbers.  I'll see if that 
> has been fixed in a library update.

Lots of arcane code and bug reporting, but I've looked at the most recent SVN 
version of gcc and see no comments in the ChangeLog regarding csqrt and a -Inf 
problem.

Before I submitted a bug report, I figured I'd have a concrete example that 
fails.  Well, for what it is worth, I tried the following little C program to 
see if it hangs:

#include <complex.h>
#include <tgmath.h>
#include <stdio.h>
int main(void) {
    _Complex double inpval = -1.0/0;
    _Complex double result = csqrt(inpval);
    printf("inpval = %lg + %lg\n", crealf(inpval), cimagf(inpval));
    printf("result = %lg + %lg\n", crealf(result), cimagf(result));
    return 0;
}

( gcc -L/usr/lib -lm csqrt_test.c ) but it doesn't hang:

address@hidden csqrt_bug]$ ./a.out
inpval = -inf + 0
result = 0 + inf
address@hidden csqrt_bug]$

So, I'm left inconclusive wondering if the compilation and library usage matches 
ultimately what Octave compilation is set up as.

Grasping at straws, but is there any possible way that Octave could be calling 
csqrt(complex value) for, say sqrt(-3), whereas it is calling sqrt(real value) 
for sqrt(-Inf)?  The latter would be a floating point exception and maybe that 
isn't handled properly?

Dan