Skip to content

Commit

Permalink
Additional speed ups of complex.pow (#23264)
Browse files Browse the repository at this point in the history
This PR speeds up complex.pow when both base and exponent are real; when
only the exponent is real; and when the base is Euler's number. These
are some pretty common cases which appear in many formulas. The speed
ups are pretty significant. According to my measurements (using the
timeit library) when both base and exponent are real the speedup is ~2x;
when only the exponent is real it is ~1.5x and when the base is Euler's
number it is ~2x.

There is no measurable difference when using other exponents which makes
sense since I refactored the code a little to reduce the total number of
branches that are needed to get to the final "fallback" branch, and
those branches have less comparisons. Anecdotally the fallback case
feels slightly faster, but the improvement is so small that I cannot
claim an improvement. If it is there it is perhaps in the order of 3 or
4%.

---------

Co-authored-by: Andreas Rumpf <rumpf_a@web.de>
  • Loading branch information
AngelEzquerra and Araq authored Jan 29, 2024
1 parent abcf45e commit 857b35c
Showing 1 changed file with 25 additions and 8 deletions.
33 changes: 25 additions & 8 deletions lib/pure/complex.nim
Original file line number Diff line number Diff line change
Expand Up @@ -249,14 +249,31 @@ func pow*[T](x, y: Complex[T]): Complex[T] =
else:
result.re = 0.0
result.im = 0.0
elif y.re == 1.0 and y.im == 0.0:
result = x
elif y.re == -1.0 and y.im == 0.0:
result = T(1.0) / x
elif y.re == 2.0 and y.im == 0.0:
result = x * x
elif y.re == 0.5 and y.im == 0.0:
result = sqrt(x)
elif y.im == 0.0:
if y.re == 1.0:
result = x
elif y.re == -1.0:
result = T(1.0) / x
elif y.re == 2.0:
result = x * x
elif y.re == 0.5:
result = sqrt(x)
elif x.im == 0.0:
# Revert to real pow when both base and exponent are real
result.re = pow(x.re, y.re)
result.im = 0.0
else:
# Special case when the exponent is real
let
rho = abs(x)
theta = arctan2(x.im, x.re)
s = pow(rho, y.re)
r = y.re * theta
result.re = s * cos(r)
result.im = s * sin(r)
elif x.im == 0.0 and x.re == E:
# Special case Euler's formula
result = exp(y)
else:
let
rho = abs(x)
Expand Down

0 comments on commit 857b35c

Please sign in to comment.