[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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

From: Richard B. Kreckel
Subject: Re: [Bug-gsl] GSL fails to honor zero's sign at branch cuts
Date: Fri, 02 Dec 2011 08:02:18 +0100
User-agent: Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv: Gecko/20111120 Icedove/3.1.16


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;

Richard B. Kreckel

reply via email to

[Prev in Thread] Current Thread [Next in Thread]