1 files changed, 44 insertions, 8 deletions
diff --git a/lib/pure/complex.nim b/lib/pure/complex.nim
index 2ea29a4f9..b48811eae 100644
--- a/lib/pure/complex.nim
+++ b/lib/pure/complex.nim
@@ -15,9 +15,6 @@
 runnableExamples:
   from std/math import almostEqual, sqrt
 
-  func almostEqual(a, b: Complex): bool =
-    almostEqual(a.re, b.re) and almostEqual(a.im, b.im)
-
   let
     z1 = complex(1.0, 2.0)
     z2 = complex(3.0, -4.0)
@@ -81,7 +78,7 @@ func abs2*[T](z: Complex[T]): T =
   ## that is the squared distance from (0, 0) to `z`.
   ## This is more efficient than `abs(z) ^ 2`.
   result = z.re * z.re + z.im * z.im
-  
+
 func sgn*[T](z: Complex[T]): Complex[T] =
   ## Returns the phase of `z` as a unit complex number,
   ## or 0 if `z` is 0.
@@ -249,10 +246,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.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)
@@ -391,6 +409,24 @@ func rect*[T](r, phi: T): Complex[T] =
   ## * `polar func<#polar,Complex[T]>`_ for the inverse operation
   complex(r * cos(phi), r * sin(phi))
 
+func almostEqual*[T: SomeFloat](x, y: Complex[T]; unitsInLastPlace: Natural = 4): bool =
+  ## Checks if two complex values are almost equal, using the
+  ## [machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon).
+  ##
+  ## Two complex values are considered almost equal if their real and imaginary
+  ## components are almost equal.
+  ##
+  ## `unitsInLastPlace` is the max number of
+  ## [units in the last place](https://en.wikipedia.org/wiki/Unit_in_the_last_place)
+  ## difference tolerated when comparing two numbers. The larger the value, the
+  ## more error is allowed. A `0` value means that two numbers must be exactly the
+  ## same to be considered equal.
+  ##
+  ## The machine epsilon has to be scaled to the magnitude of the values used
+  ## and multiplied by the desired precision in ULPs unless the difference is
+  ## subnormal.
+  almostEqual(x.re, y.re, unitsInLastPlace = unitsInLastPlace) and
+  almostEqual(x.im, y.im, unitsInLastPlace = unitsInLastPlace)
 
 func `$`*(z: Complex): string =
   ## Returns `z`'s string representation as `"(re, im)"`.