From bug-octave-request at bevo dot che dot wisc dot edu Tue Jan 22 16:11:55 2002 Subject: Re: IMPORTANT: Bug in "arg" function From: Paul Kienzle To: robher at adinet dot com dot uy Cc: bug-octave at bevo dot che dot wisc dot edu Date: Tue, 22 Jan 2002 17:11:51 -0500 Fundamentally, arg(0) is ill-defined regardless of what the sign is on the 0. Close to 0, the phase calculation is highly sensitive to error, so you should be setting the phase of all values with really small magnitude by hand. According to http://www.lysator.liu.se/c/rat/d5.html, atan2 is not standardized when both arguments are 0: The atan2 function is modelled after FORTRAN's. It is described in terms of arctany/x for simplicity; the Committee did not wish to complicate the descriptions by specifying in detail how the determine the appropriate quadrant, since that should be obvious from normal mathematical convention. atan2(y,x) is well-defined and finite, even when x is 0; the one ambiguity occurs when both arguments are 0, because at that point any value in the range of the function could logically be selected. Since valid reasons can be advanced for all the different choices that have been in this situation by various implements, the Standard preserves the implementor's freedom to return an arbitrary well-defined value such as 0, to report a domain error, or to return an IEEE NaN code. According to IEEE 754/854, 1/-Inf -> -0 and 1/+Inf -> +0, so in the limit you have -0 == +0. If you are using Inf as an infinity, then you should have -0 == +0. If you are using Inf as some number bigger than REAL_MAX, then you should have -0 < +0. Both situations occur in practice: 1/0 is infinite, 10/REAL_MIN is bigger than REAL_MAX but still finite. In the second situation, you really want atan2(0,0)=0 and atan2(0,-0)=pi. In the first, you want them both to be NaN. Since differentiating the results is more flexible, I agree with glibc. Too bad other libc's (e.g., IRIX) return 0 for both +0 and -0. Should octave set a standard where others agree not to have one? Paul Kienzle pkienzle at users dot sf dot net On Tue, Jan 22, 2002 at 05:30:08PM -0300, Roberto Hernandez wrote: > Paul Kienzle wrote: > > > After playing around in the debugger for a while, it seems that the > problem > > is that (-0 < 0) is not true. > > ***** Awesome, you figured it out! *************** > > > Looking at the value X(3), the debugger says that this is -0+0i. When > > evaluated as a vector, the function complex::arg in the standard > > C++ library returns atan2 (imag (x), real (x)), which apparently > recognizes > > -0 as a negative number. However, when evaluated as a scalar, Octave > > automatically converts the type to real(X(3)) since the complex portion > > is 0. This is then evaluated in liboctave/lo-mappers.cc using the > > following: > > > > double > > arg (double x) > > { > > if (x < 0.0) > > return M_PI; > > else > > #if defined (HAVE_ISNAN) > > return xisnan (x) ? octave_NaN : 0.0; > > #else > > return 0.0; > > #endif > > } > > > > Since x < 0.0 is false for x=-0.0, this evaluates to 0.0 rather than > > M_PI. You get consistent behaviour if you replace > > if (x < 0.0) > > with > > double y = copysign(1.0,x); > > if (y < 0.0) > > > OK, I'm a bit confused. In the example originally posted, the correct > answer was the one displayed by evaluating each zero element in the > vector as a scalar. In other words, in this case it's good that the > function "arg" evaluates (-0.0 < 0.0) as false. Shouldn't the behavior > be made consistent the other way around? > > I guess my confusion lies in not understanding what -0.0 is supposed to > mean. The way I originally noticed the "arg" bug was while doing some > FFT processing. The transform of a simple sinewave, which is supposed to > have just two peaks in the phase plot looked like a sawtooth. > > After some digging I figured out that the FFT was returning some > residual values where it was supposed to return a 0. In this case it was > OK that the "arg" function returned weird phase values, because the > input vector was messed up. > > So my next action was to write a function to round all the "small" > values to zero, that is to _ABSOLUTE_ zero. Then I passed the "clean" > vector to "arg" and that's when I noticed the weird behavior. In this > case it shouldn't be OK that "arg" returns PI for some zero values, > because they are absolute zero and not some infinitesimal value. > > In summary, when does the -0.0 come into play? Is it supposed to > represent an infinitesimal negative value? Does (-0.0 == 0.0) evaluate > as true? > > IMHO we should adopt some criteria as to what is the "CORRECT" thing for > "arg" to do when passed a -0.0 value, whether real or imaginary. It > occurs to me that if the function complex::arg is passed 0-0i it > would return -PI/2, which is incorrect. > > Opinions? > > > Thanks again, > > Roberto > > ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------