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

[Charles R Harris]

On Sat, Jul 19, 2008 at 7:13 AM, Pauli Virtanen <pav at 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.
>
>
Agreed.


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

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20080719/16219ff8/attachment.html>