[Numpy-discussion] Branch cuts, inf, nan, C99 compliance

Charles R Harris charlesr.harris@gmail....
Sat Jul 19 12:05:19 CDT 2008

On Sat, Jul 19, 2008 at 7:13 AM, Pauli Virtanen <pav@iki.fi> wrote:

> Hi all,
> Re: Ticket 854.
> I wrote tests for the branch cuts for all complex arc* functions
> in umathmodule. It turns out that all except arccosh were OK.
> The formula for arcsinh was written in a non-standard form with
> an unnecessary nc_neg, but this didn't affect the results.
> I also wrote tests for checking values of the functions at infs and nans.
> A patch for all of this is attached, with all currently non-passing
> tests marked as skipped. I'd like to commit this if there are no
> objections.
> Another thing I noticed is that the present implementations of
> the complex functions are naive, so they over- and underflow earlier
> than necessary:
> >>> np.arcsinh(1e8)
> 19.113827924512311
> >>> np.arcsinh(1e8 + 0j)
> (inf-0j)
> >>> np.arcsinh(-1e8 + 0j)
> (-19.113827924512311-0j)
> This particular thing in arcsinh occurs because of loss of precision
> in intermediate stages. (In the version in my patch this loss of precision
> is still present.)
> It would be nice to polish these up. BTW, are there obstacles to using
> the C99 complex functions when they are available? This would avoid
> quite a bit of drudgework... As an alternative, we can probably steal
> better implementations from Python's recently polished cmathmodule.c

The main problem is determining when they are available and if they cover
all the needed precisions. Since we will need standalone implementations on
some platforms anyway, I am inclined towards stealing from cmathmodule.c if
it offers improved code for some of the functions.

>    ***
> Then comes a descent into pedantry: Numpy complex functions are not
> C99 compliant in handling of the signed zero, inf, and nan. I don't
> know whether we should comply with C99, but it would be nicer to have
> these handled in a controlled way rather than as a byproduct of the
> implementation chosen.

> 1)
> The branch cuts for sqrt and arc* don't respect the negative zero:
> >>> a = 1 + 0j
> >>> np.sqrt(-a)
> 1j
> >>> np.sqrt(-1 + 0j)
> 1j
> The branch cut of the logarithm however does:
> >>> np.log(-a)
> -3.1415926535897931j
> >>> np.log(-1 + 0j)
> 3.1415926535897931j
> All complex functions in the C99 standard respect the negative zero
> in their branch cuts. Do we want to follow?

Hmm. I think so, to the extent we can. This might lead to some unexpected
results, but that is what happens for arguments near the branch cuts. Do it
and document it.

> I don't know how to check what Octave and Matlab do regarding this,
> since I haven't figured out how to place a negative zero in complex
> numbers in these languages. But at least in practice these languages
> appear not to respect the sign of zero.
> > a = 1 + 0j
> > log(-a)
> ans =  0.000000000000000 + 3.141592653589793i
> > log(-1)
> ans =  0.000000000000000 + 3.141592653589793i
> 2)
> The numpy functions in general don't return C99 compliant results
> at inf or nan. I wrote up some tests for checking these.
> Do we want to fix these?

I'd say yes. That way we can refer to the C99 standard to document numpy

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20080719/16219ff8/attachment.html 

More information about the Numpy-discussion mailing list