Faster, more correct complex^complex

[stevengj]

This PR fixes the incorrect behaviors for complex^complex identified in #24515, and it also makes the code significantly faster without (as far as I can tell) sacrificing accuracy. On my machine, it is around 60% faster for complex^complex and 120% faster for real^complex in double precision.

The old code had two completely separate implementations, one for floating-point types and one for other types, despite the fact that both produced floating-point results. The floating-point version was based on z^p = exp(p * log(z)), whereas the other version first converted z to polar form and then exponentiated. The latter approach seems to be significantly faster and no less accurate, so I now use that in all cases (with various special-case optimizations for real z and/or p). By unifying the implementations, the code is also significantly shorter.

Marked as breaking only because this throws exceptions in some cases where the old code would have silently returned NaNs. Returns NaN as before.

See also the discussions in #2891, #3246, 06530b6.

To do:

• Tests: need to cover all of the branches and edge cases.
• Fixes/tests for Inf and NaN cases
• Decide whether to throw DomainError for (0.0+0.0im)^-1.0 and other cases ala NaN vs wild (or, what's a DomainError, really?) #5234

[oscardssmith]

Any chance of getting nice error plots to make sure that error stays below x ulps?

[StefanKarpinski]

Maybe write it as exp2(31) instead? It is jarring to see integer types in here.

[stevengj]

But typemax(Int32) is literally what we want to prevent overflow on 32-bit machines, and to have consistent behavior on 64-bit. The whole point of this if statement is to convert to an integer type.

[stevengj]

@oscardssmith, some plots of the old and new relative errors for z^w and w^z for a variety of z (in double precision) are shown below: the accuracy seems very similar. I should also note that the polar formula (used in my new code) is how CPython computes it and is also how R computes it, so we are in pretty good company.

The main trick is to get the edge cases (especially the sign of zero imaginary parts) correct, especially when you want to take advantage of real z and/or p for z^p.

[stevengj]

@StefanKarpinski, is the policy these days to throw DomainError rather than constructing NaN explicitly, if it is practical to do so? inv(0.0+0.0im) (and hence (0.0+0.0im)^-1) doesn't throw an error, but I guess that is because it would slow it down too much to check for this?

I'm a little confused because it doesn't seem like #5234 was ever resolved.

[giordano]

What are the white points in the plots?

[stevengj]

@giordano, the white points are where the error is "zero" (i.e. the result, by luck, is exactly rounded), so log of the error is –∞.

Note that the error is unlikely to be actually zero for the white points, it is just < 1ulp, so probably they should be bluish or greenish. Here are the plots re-done with the true error (i.e. I compute the error and then round from BigFloat to Float64 rather than the other way around), making sure the color scales match:

[giordano]

Uhm, to me the current implementation seems to have more white points, on the other hand the proposed one looks more blueish, which is good anyway.

[stevengj]

If I change this to throw exceptions rather than creating NaNs from non-NaN inputs, that will technically be a breaking change. Again, I'm not sure what the desired behavior is these days?

[stevengj]

By the way, one argument against throwing DomainError from this function: if we want to catch all cases where non-NaN inputs can generate NaN results, there are a whole bunch of additional overflow circumstances etc. that need to be checked.

Also, z^p for p::Integer does not throw, because power_by_squaring does no NaN checks. e.g. (Inf+Inf*im)^3 produces NaN + NaN*im.

Are we okay with throwing DomainError in some NaN-producing cases but not all? The alternatives are:

• add a lot more checks (making the code more messy and possibly slower) … it's pretty hard to be sure we've covered everything, too

• or never throw (just silently return NaNs).

[stevengj]

After sleeping on it, I've decided to revert the breaking change of throwing DomainError. Rationale:

• This makes the PR strictly more conservative, so throwing errors can always be added in a separate PR.

• It didn't make sense to me to throw DomainError sometimes and generate NaN sometimes. It should be one or the other.

• Finding all cases where non-NaN inputs produce NaN would greatly complicate the code, and would make it inconsistent with integer powers (power_by_squaring).

[JeffBezanson]

We don't have a general rule that getting NaNs from non-NaN inputs should be a DomainError. That pattern was introduced just for a couple functions where that check happens to correctly identify DomainError cases.

[stevengj]

CI looks good except AppVeyor, where the error Can not open the file as [PE] archive looks unrelated.

[pkofod]

(i.e. I compute the error and then round from BigFloat to Float64 rather than the other way around)

Just to be sure does this mean that you are converting the Float64 result to BigFloat, subtract a reference BigFloat result, divide by (?), and convert back to Float64?

[stevengj]

Just to be sure does this mean that you are converting the Float64 result to BigFloat, subtract a reference BigFloat result, divide by (?), and convert back to Float64?

Essentially yes. To compute the relative error for the Complex{Float64} value x = z^p in the updated plots, I first get the "exact" result e = big(z)^big(p), then compute Float64(abs(x - e) / abs(e)). The x-e I think may call a BigFloat-Float64 MPFR function to do the subtraction directly without first converting x to Complex{BigFloat}, however.

[stevengj]

Should be good to merge?

[KristofferC]

@nanosoldier runbenchmarks(ALL, vs = ":master")

[stevengj]

Is nanosoldier working? It's been "running on node 3" for a day now.

[KristofferC]

The only thing we can conclude is that this PR brings some serious performance regressions ;)

[stevengj]

Bump. Seems like there is no point in running nanosoldier right now...

[ararslan]

@nanosoldier runbenchmarks(ALL, vs=":master")

[nanosoldier]

Something went wrong when running your job:

NanosoldierError: failed to run benchmarks against primary commit: failed process: Process(`sudo cset shield -e su nanosoldier -- -c ./benchscript.sh`, ProcessExited(1)) [1]

Logs and partial data can be found here
cc @ararslan

[stevengj]

Looks like nanosoldier is not working; seems like an unrelated ProcessExited(1) exception since it occurs before of the be benchmarks start running.

[ararslan]

Okay, not sure yet what the deal is. In the meantime you can benchmark manually with BenchmarkTools, since the macros work on 0.7 again now.

[stevengj]

Let's see if nanosoldier works now:

@nanosoldier runbenchmarks(ALL, vs=":master")

(Though it has been so noisy recently, even when it works, that the results have been hard to make use of.)

[stevengj]

@ararslan, what's up with nanosoldier?

[ararslan]

It's been acting up a bit lately. I think in this case I restarted the server after you had triggered a run, which means that the status for that run didn't get updated. I'll try retriggering.

@nanosoldier runbenchmarks(ALL, vs=":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @ararslan

[stevengj]

The nanosoldier regressions look like noise.

(It doesn't look like there are any benchmarks in BaseBenchmarks that exercise complex powers.)

[stevengj]

Seems ready to merge?

[stevengj]

Yay, CI is green.