Re: [Bug-gsl] GSL fails to honor zero's sign at branch cuts

[Richard B. Kreckel]

Hi!

On 11/15/2011 11:23 PM, Brian Gough wrote:
> At Sun, 13 Nov 2011 00:23:32 +0100,
> Richard B. Kreckel wrote:
> > With that patch, some checks fail. Where they fail it should be
> > sufficient to flip the sign of the zero real or imaginary part in the
> > argument. (Of course, if one wishes, one can add new test vectors for
> > the opposite signed zero.)
>
> Hello Richard.  Thanks for the patch.  I'm happy to apply this and
> follow the general signed zero convention

Thanks. For the sake of completeness, let me mention that the current 
GSL already treats the branch cut of the complex logarithm correctly: It 
distinguishes between log(x+0i) and log(x-0i) by virtue of the relative 
sign of x and y in atan2(x,y). Though, a test involving -0.0 appears to 
be missing from complex/results.h.

> -- I'm not in a position to
> check it myself so I'll just go with whatever you submit, but I just
> wanted to confirm some things:
>
> - in the manual we say "for multiple-valued functions the branch cuts
> have been chosen to follow the conventions of Abramowitz and Stegun in
> the Handbook of Mathematical Functions. The functions return
> principal values which are the same as those in GNU Calc, which in
> turn are the same as those in Common Lisp, The Language (Second
> Edition) and the HP-28/48 series of calculators."

The *positions* of the branch cuts are the same in CLTL2, C99, 
Abromowitz and Stegun, etc.

However, a system with signed zero cab give different results as a 
system without signed zero *on* the branch cut! For example, take the 
branch cut of asin starting at 1 and running along the positive real 
axis towards infinity. A system with signed zero can distinguish between 
asin(5+0i) and asin(5-0i) and pick the appropriate branch cut continuous 
with quadrant I or quadrant IV, respectively. A system without signed 
zero cannot distinguish and, according to the choses the branch cut such 
that it is continuous with quadrant IV (consistent with the rule of 
CCC). So, for asin(5+0i) it will return the same value as a system with 
signed zero would for asin(5-0i).

> Does this still hold?  I generated the test values with GNU Calc
> (complex/test.el) so there's some kind of inconsistency if the tests
> are failing (maybe because GNU Calc does not have signed zero?).
> Maybe they should be generated with a common lisp implementation with
> signed zero instead.

Right. I propose to generate them using the MPC library 
<mpc.multiprecision.org>. The MPC developers did a fine job in taking 
care that the branch cuts are the same as in C99. And I just checked 
MPC's log, sqrt, asin[h], acos[h], and atan[h] and can confirm that they 
comply.

> - we aim to support old C89 compilers.  Is there a portable way to do
> a signbit() function in C89?

Not without making assumptions about endianness of floating-point 
variables, I'm afraid. But caring about this it is straightfoward. I 
assume we cannot rely on the <endian.h> header either. Does GSL already 
have configure scripts that test the endianness of double and how it 
maps to ints? It would then have to be done similar to this:

int signbit(double x)
{
    union { double d; int i[2]; } u = { d: x };
    return u.i[1] < 0;
}

Bye!
  -richy.
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>