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

[Richard B. Kreckel]

Systems that support signed zero should be careful about the value on 
the branch cuts ("slits"). To quote William Kahan [0]:

"Since all the slits in question lie on either the real or the imaginary 
axis, every point z on a slit is represented in at least two ways, at 
least once with a +0 and at least once with a -0 for whichever of the 
real and imaginary parts of z vanishes. Benignly, ambiguity in z at a 
discontinuity of f(z) permits f(z) to be defined formally continuously, 
except possibly at the ends of some slits, by continuation from inside 
the principal domain. This continuity goes beyond mere formality. By 
analytic continuation, the domain of each of our nine elementary 
functions f(z) extends until it fills out a Riemann Surface; think of 
this surface as a multiple covering wrapped like a bandage around the 
Riemann Sphere and mapped onto it continuously by f. To construct f's 
principal domain, cut the bandage along the slit(s) and discard all but 
one layer covering the sphere. That layer is a closed surface mapped by 
f continuously onto a subset of the sphere. The shadow of that layer 
projected down upon the sphere is the principal domain; it consists of 
the whole sphere, but with slit(s) covered twice. That is why we wish to 
represent slits ambiguously."

The basic convention (adopted in C99, CLTL2, etc. and implemented in 
GCC, GLibC, CLisp, etc.) for branch cuts with only one finite endpoint 
is counter-clockwise continuity (CCC): "the function is continuous as 
the cut is approached coming around the finite endpoint of the cut in a 
counter clock wise direction" (wording from [1]).

GSL falls into the category of systems that should do this, since it 
uses IEEE 754 double precision for constructing complex numbers. 
However, GSL does not distinguish between +0 and -0 on the branch cuts!

I propose to fix that for two reasons:
1) Distinguishing between +0 and -0 is a superior concept.
2) Fixing this harmonizes GSL with other systems that support both 
signed zero and complex numbers.

I attach a patch into that direction. It fixes sqrt, arcsin, arccos, 
arctan, arcsinh, arccosh, and arctanh.

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.)

Thanks for applying.

  -richy.

[0] William Kahan: "Branch Cuts for Complex Elementary Functions, or 
Much Ado About Nothing's Sign Bit"
[1] ISO/IEC 9899:1999 (E): "Programming Languages: C"
--
Richard B. Kreckel
<http://www.ginac.de/~kreckel/>

gsl.patch
Description: Text Data