Implements Kahan's versions of elementary complex functions

[jiahao]

I've implemented Kahan's versions of elementary complex functions, paying particular attention to behavior near branch cuts.

Why ?

The behavior of complex elementary functions has come up in recent issues (#2889 #2843 #2845 #2875) and it was probably high time to have implementations of these functions that paid proper attention to branch cut issues.
Hopefully this is a good starting point to get the complex arithmetic done up properly in Julia.

What

This commit provides the functions described in the Reference, implementing these algorithms verbatim (except where noted).

Ref: Branch cuts for complex elementary functions

• List of functions affected: sqrt, ^, log, asin, acos, atan, tanh, asinh, acosh, atanh
• Includes also ssqs to compute the mod-squared of a complex number scaled by a power of 2 to avoid overflow/underflow. ssqs is not exported.
• the part about dealing with unsigned zeros in atanh in the reference is not implemented, this is not necessary with standards-compliant floating point.

Issues

• The ^ function fails tests so is disabled for now.
• The log function works only for Complex{FloatingPoint}. Explicitly casting {T<:Integer}Complex{T} to floating-point Complex causes tests to fail.
• Implementations need to be checked for correctness
• Behaviors need to be checked
• All tests in Fix complex function branch cut behaviour #2875 need to pass
• Additional testing especially near branch cuts

@alanedelman: "Nobody pays attention to branch cuts even for scalar complex functions, let alone matrix-valued functions"

cc: @ViralBShah @JeffBezanson @simonbyrne @bsxfan @andreasnoackjensen

[JeffBezanson]

Are all the tests in complex.jl updated to reflect our intended answers?

[jiahao]

Not yet; however, those are not the tests that are failing atm

[JeffBezanson]

I believe sqrt and log were already written to avoid overflow; can you provide cases where they got a worse answer than the new versions?

[jiahao]

Not offhand but I am running an exhaustive search along the branch cut areas to see. Thus far I've discovered some bugs in my implementations.

log is the "easy" case and sqrt is pretty well studied so I wouldn't be surprised if they were already optimal in a sense. The inverse functions are tricky because they have funny branch cuts; Kahan discusses them at length and cautions against how the "obvious" formulae don't work.

[simonbyrne]

One example of an overflow problem in the current version:

julia> tanh(1e3)
1.0
julia> tanh(complex(1e3,0.0))
NaN + 0.0im

Interestingly, numpy has a similar problem

I haven't tried it, but this patch would seem to fix it.

[jiahao]

I can confirm that tanh does not overflow in Kahan's version.

[ViralBShah]

Bump. It would be nice to get this merged.

[jiahao]

Unfortunately I have to wind up my other MIT commitments in short order and the earliest I will be able to look at this more seriously will be after June 15.

[ViralBShah]

Good to know that we will have more of your time in a few weeks!

[JeffBezanson]

The new versions of asin, acos, atan, tanh, asinh, acosh, atanh are pretty clearly improvements so I'd be ok merging those.

[jiahao]

I'm fine with letting someone take over this if he/she felt like it.

[JeffBezanson]

I have test-complex passing on this branch now. @simonbyrne what do you think?

[simonbyrne]

Fantastic work: based on my limited knowledge, it looks pretty good. If the tests are passing, you can add complex to the list in runtests.jl.