diff options
Diffstat (limited to 'lib/pure/complex.nim')
| -rw-r--r-- | lib/pure/complex.nim | 107 |
1 files changed, 87 insertions, 20 deletions
diff --git a/lib/pure/complex.nim b/lib/pure/complex.nim index b9371c1e1..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) @@ -36,7 +33,7 @@ runnableExamples: {.push checks: off, line_dir: off, stack_trace: off, debugger: off.} # the user does not want to trace a part of the standard library! -import std/math +import std/[math, strformat] type Complex*[T: SomeFloat] = object @@ -82,6 +79,13 @@ func abs2*[T](z: Complex[T]): T = ## 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. + let a = abs(z) + if a != 0: + result = z / a + func conjugate*[T](z: Complex[T]): Complex[T] = ## Returns the complex conjugate of `z` (`complex(z.re, -z.im)`). result.re = z.re @@ -156,18 +160,7 @@ func `/`*[T](x: T; y: Complex[T]): Complex[T] = func `/`*[T](x, y: Complex[T]): Complex[T] = ## Divides two complex numbers. - var r, den: T - if abs(y.re) < abs(y.im): - r = y.re / y.im - den = y.im + r * y.re - result.re = (x.re * r + x.im) / den - result.im = (x.im * r - x.re) / den - else: - r = y.im / y.re - den = y.re + r * y.im - result.re = (x.re + r * x.im) / den - result.im = (x.im - r * x.re) / den - + x * conjugate(y) / abs2(y) func `+=`*[T](x: var Complex[T]; y: Complex[T]) = ## Adds `y` to `x`. @@ -253,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) @@ -395,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)"`. @@ -403,4 +435,39 @@ func `$`*(z: Complex): string = result = "(" & $z.re & ", " & $z.im & ")" +proc formatValueAsTuple(result: var string; value: Complex; specifier: string) = + ## Format implementation for `Complex` representing the value as a (real, imaginary) tuple. + result.add "(" + formatValue(result, value.re, specifier) + result.add ", " + formatValue(result, value.im, specifier) + result.add ")" + +proc formatValueAsComplexNumber(result: var string; value: Complex; specifier: string) = + ## Format implementation for `Complex` representing the value as a (RE+IMj) number + ## By default, the real and imaginary parts are formatted using the general ('g') format + let specifier = if specifier.contains({'e', 'E', 'f', 'F', 'g', 'G'}): + specifier.replace("j") + else: + specifier.replace('j', 'g') + result.add "(" + formatValue(result, value.re, specifier) + if value.im >= 0 and not specifier.contains({'+', '-'}): + result.add "+" + formatValue(result, value.im, specifier) + result.add "j)" + +proc formatValue*(result: var string; value: Complex; specifier: string) = + ## Standard format implementation for `Complex`. It makes little + ## sense to call this directly, but it is required to exist + ## by the `&` macro. + ## For complex numbers, we add a specific 'j' specifier, which formats + ## the value as (A+Bj) like in mathematics. + if specifier.len == 0: + result.add $value + elif 'j' in specifier: + formatValueAsComplexNumber(result, value, specifier) + else: + formatValueAsTuple(result, value, specifier) + {.pop.} |