mirror of
https://github.com/nim-lang/Nim.git
synced 2025-12-28 17:04:41 +00:00
d8e1504ed1
This adds a version of `almostEqual` (which was already available for floats) thata works with `Complex[SomeFloat]`. Proof that this is needed is that the first thing that the complex.nim runnable examples block did before this commit was define (an incomplete) `almostEqual` function that worked with complex values.
474 lines
14 KiB
Nim
474 lines
14 KiB
Nim
#
|
|
#
|
|
# Nim's Runtime Library
|
|
# (c) Copyright 2010 Andreas Rumpf
|
|
#
|
|
# See the file "copying.txt", included in this
|
|
# distribution, for details about the copyright.
|
|
#
|
|
|
|
## This module implements complex numbers
|
|
## and basic mathematical operations on them.
|
|
##
|
|
## Complex numbers are currently generic over 64-bit or 32-bit floats.
|
|
|
|
runnableExamples:
|
|
from std/math import almostEqual, sqrt
|
|
|
|
let
|
|
z1 = complex(1.0, 2.0)
|
|
z2 = complex(3.0, -4.0)
|
|
|
|
assert almostEqual(z1 + z2, complex(4.0, -2.0))
|
|
assert almostEqual(z1 - z2, complex(-2.0, 6.0))
|
|
assert almostEqual(z1 * z2, complex(11.0, 2.0))
|
|
assert almostEqual(z1 / z2, complex(-0.2, 0.4))
|
|
|
|
assert almostEqual(abs(z1), sqrt(5.0))
|
|
assert almostEqual(conjugate(z1), complex(1.0, -2.0))
|
|
|
|
let (r, phi) = z1.polar
|
|
assert almostEqual(rect(r, phi), z1)
|
|
|
|
{.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, strformat]
|
|
|
|
type
|
|
Complex*[T: SomeFloat] = object
|
|
## A complex number, consisting of a real and an imaginary part.
|
|
re*, im*: T
|
|
Complex64* = Complex[float64]
|
|
## Alias for a complex number using 64-bit floats.
|
|
Complex32* = Complex[float32]
|
|
## Alias for a complex number using 32-bit floats.
|
|
|
|
func complex*[T: SomeFloat](re: T; im: T = 0.0): Complex[T] =
|
|
## Returns a `Complex[T]` with real part `re` and imaginary part `im`.
|
|
result.re = re
|
|
result.im = im
|
|
|
|
func complex32*(re: float32; im: float32 = 0.0): Complex32 =
|
|
## Returns a `Complex32` with real part `re` and imaginary part `im`.
|
|
result.re = re
|
|
result.im = im
|
|
|
|
func complex64*(re: float64; im: float64 = 0.0): Complex64 =
|
|
## Returns a `Complex64` with real part `re` and imaginary part `im`.
|
|
result.re = re
|
|
result.im = im
|
|
|
|
template im*(arg: typedesc[float32]): Complex32 = complex32(0, 1)
|
|
## Returns the imaginary unit (`complex32(0, 1)`).
|
|
template im*(arg: typedesc[float64]): Complex64 = complex64(0, 1)
|
|
## Returns the imaginary unit (`complex64(0, 1)`).
|
|
template im*(arg: float32): Complex32 = complex32(0, arg)
|
|
## Returns `arg` as an imaginary number (`complex32(0, arg)`).
|
|
template im*(arg: float64): Complex64 = complex64(0, arg)
|
|
## Returns `arg` as an imaginary number (`complex64(0, arg)`).
|
|
|
|
func abs*[T](z: Complex[T]): T =
|
|
## Returns the absolute value of `z`,
|
|
## that is the distance from (0, 0) to `z`.
|
|
result = hypot(z.re, z.im)
|
|
|
|
func abs2*[T](z: Complex[T]): T =
|
|
## Returns the squared absolute value of `z`,
|
|
## 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.
|
|
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
|
|
result.im = -z.im
|
|
|
|
func inv*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the multiplicative inverse of `z` (`1/z`).
|
|
conjugate(z) / abs2(z)
|
|
|
|
func `==`*[T](x, y: Complex[T]): bool =
|
|
## Compares two complex numbers for equality.
|
|
result = x.re == y.re and x.im == y.im
|
|
|
|
func `+`*[T](x: T; y: Complex[T]): Complex[T] =
|
|
## Adds a real number to a complex number.
|
|
result.re = x + y.re
|
|
result.im = y.im
|
|
|
|
func `+`*[T](x: Complex[T]; y: T): Complex[T] =
|
|
## Adds a complex number to a real number.
|
|
result.re = x.re + y
|
|
result.im = x.im
|
|
|
|
func `+`*[T](x, y: Complex[T]): Complex[T] =
|
|
## Adds two complex numbers.
|
|
result.re = x.re + y.re
|
|
result.im = x.im + y.im
|
|
|
|
func `-`*[T](z: Complex[T]): Complex[T] =
|
|
## Unary minus for complex numbers.
|
|
result.re = -z.re
|
|
result.im = -z.im
|
|
|
|
func `-`*[T](x: T; y: Complex[T]): Complex[T] =
|
|
## Subtracts a complex number from a real number.
|
|
result.re = x - y.re
|
|
result.im = -y.im
|
|
|
|
func `-`*[T](x: Complex[T]; y: T): Complex[T] =
|
|
## Subtracts a real number from a complex number.
|
|
result.re = x.re - y
|
|
result.im = x.im
|
|
|
|
func `-`*[T](x, y: Complex[T]): Complex[T] =
|
|
## Subtracts two complex numbers.
|
|
result.re = x.re - y.re
|
|
result.im = x.im - y.im
|
|
|
|
func `*`*[T](x: T; y: Complex[T]): Complex[T] =
|
|
## Multiplies a real number with a complex number.
|
|
result.re = x * y.re
|
|
result.im = x * y.im
|
|
|
|
func `*`*[T](x: Complex[T]; y: T): Complex[T] =
|
|
## Multiplies a complex number with a real number.
|
|
result.re = x.re * y
|
|
result.im = x.im * y
|
|
|
|
func `*`*[T](x, y: Complex[T]): Complex[T] =
|
|
## Multiplies two complex numbers.
|
|
result.re = x.re * y.re - x.im * y.im
|
|
result.im = x.im * y.re + x.re * y.im
|
|
|
|
func `/`*[T](x: Complex[T]; y: T): Complex[T] =
|
|
## Divides a complex number by a real number.
|
|
result.re = x.re / y
|
|
result.im = x.im / y
|
|
|
|
func `/`*[T](x: T; y: Complex[T]): Complex[T] =
|
|
## Divides a real number by a complex number.
|
|
result = x * inv(y)
|
|
|
|
func `/`*[T](x, y: Complex[T]): Complex[T] =
|
|
## Divides two complex numbers.
|
|
x * conjugate(y) / abs2(y)
|
|
|
|
func `+=`*[T](x: var Complex[T]; y: Complex[T]) =
|
|
## Adds `y` to `x`.
|
|
x.re += y.re
|
|
x.im += y.im
|
|
|
|
func `-=`*[T](x: var Complex[T]; y: Complex[T]) =
|
|
## Subtracts `y` from `x`.
|
|
x.re -= y.re
|
|
x.im -= y.im
|
|
|
|
func `*=`*[T](x: var Complex[T]; y: Complex[T]) =
|
|
## Multiplies `x` by `y`.
|
|
let im = x.im * y.re + x.re * y.im
|
|
x.re = x.re * y.re - x.im * y.im
|
|
x.im = im
|
|
|
|
func `/=`*[T](x: var Complex[T]; y: Complex[T]) =
|
|
## Divides `x` by `y` in place.
|
|
x = x / y
|
|
|
|
|
|
func sqrt*[T](z: Complex[T]): Complex[T] =
|
|
## Computes the
|
|
## ([principal](https://en.wikipedia.org/wiki/Square_root#Principal_square_root_of_a_complex_number))
|
|
## square root of a complex number `z`.
|
|
var x, y, w, r: T
|
|
|
|
if z.re == 0.0 and z.im == 0.0:
|
|
result = z
|
|
else:
|
|
x = abs(z.re)
|
|
y = abs(z.im)
|
|
if x >= y:
|
|
r = y / x
|
|
w = sqrt(x) * sqrt(0.5 * (1.0 + sqrt(1.0 + r * r)))
|
|
else:
|
|
r = x / y
|
|
w = sqrt(y) * sqrt(0.5 * (r + sqrt(1.0 + r * r)))
|
|
|
|
if z.re >= 0.0:
|
|
result.re = w
|
|
result.im = z.im / (w * 2.0)
|
|
else:
|
|
result.im = if z.im >= 0.0: w else: -w
|
|
result.re = z.im / (result.im + result.im)
|
|
|
|
func exp*[T](z: Complex[T]): Complex[T] =
|
|
## Computes the exponential function (`e^z`).
|
|
let
|
|
rho = exp(z.re)
|
|
theta = z.im
|
|
result.re = rho * cos(theta)
|
|
result.im = rho * sin(theta)
|
|
|
|
func ln*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the
|
|
## ([principal value](https://en.wikipedia.org/wiki/Complex_logarithm#Principal_value)
|
|
## of the) natural logarithm of `z`.
|
|
result.re = ln(abs(z))
|
|
result.im = arctan2(z.im, z.re)
|
|
|
|
func log10*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the logarithm base 10 of `z`.
|
|
##
|
|
## **See also:**
|
|
## * `ln func<#ln,Complex[T]>`_
|
|
result = ln(z) / ln(10.0)
|
|
|
|
func log2*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the logarithm base 2 of `z`.
|
|
##
|
|
## **See also:**
|
|
## * `ln func<#ln,Complex[T]>`_
|
|
result = ln(z) / ln(2.0)
|
|
|
|
func pow*[T](x, y: Complex[T]): Complex[T] =
|
|
## `x` raised to the power of `y`.
|
|
if x.re == 0.0 and x.im == 0.0:
|
|
if y.re == 0.0 and y.im == 0.0:
|
|
result.re = 1.0
|
|
result.im = 0.0
|
|
else:
|
|
result.re = 0.0
|
|
result.im = 0.0
|
|
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)
|
|
theta = arctan2(x.im, x.re)
|
|
s = pow(rho, y.re) * exp(-y.im * theta)
|
|
r = y.re * theta + y.im * ln(rho)
|
|
result.re = s * cos(r)
|
|
result.im = s * sin(r)
|
|
|
|
func pow*[T](x: Complex[T]; y: T): Complex[T] =
|
|
## The complex number `x` raised to the power of the real number `y`.
|
|
pow(x, complex[T](y))
|
|
|
|
|
|
func sin*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the sine of `z`.
|
|
result.re = sin(z.re) * cosh(z.im)
|
|
result.im = cos(z.re) * sinh(z.im)
|
|
|
|
func arcsin*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse sine of `z`.
|
|
result = -im(T) * ln(im(T) * z + sqrt(T(1.0) - z*z))
|
|
|
|
func cos*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the cosine of `z`.
|
|
result.re = cos(z.re) * cosh(z.im)
|
|
result.im = -sin(z.re) * sinh(z.im)
|
|
|
|
func arccos*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse cosine of `z`.
|
|
result = -im(T) * ln(z + sqrt(z*z - T(1.0)))
|
|
|
|
func tan*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the tangent of `z`.
|
|
result = sin(z) / cos(z)
|
|
|
|
func arctan*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse tangent of `z`.
|
|
result = T(0.5)*im(T) * (ln(T(1.0) - im(T)*z) - ln(T(1.0) + im(T)*z))
|
|
|
|
func cot*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the cotangent of `z`.
|
|
result = cos(z)/sin(z)
|
|
|
|
func arccot*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse cotangent of `z`.
|
|
result = T(0.5)*im(T) * (ln(T(1.0) - im(T)/z) - ln(T(1.0) + im(T)/z))
|
|
|
|
func sec*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the secant of `z`.
|
|
result = T(1.0) / cos(z)
|
|
|
|
func arcsec*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse secant of `z`.
|
|
result = -im(T) * ln(im(T) * sqrt(1.0 - 1.0/(z*z)) + T(1.0)/z)
|
|
|
|
func csc*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the cosecant of `z`.
|
|
result = T(1.0) / sin(z)
|
|
|
|
func arccsc*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse cosecant of `z`.
|
|
result = -im(T) * ln(sqrt(T(1.0) - T(1.0)/(z*z)) + im(T)/z)
|
|
|
|
func sinh*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the hyperbolic sine of `z`.
|
|
result = T(0.5) * (exp(z) - exp(-z))
|
|
|
|
func arcsinh*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse hyperbolic sine of `z`.
|
|
result = ln(z + sqrt(z*z + 1.0))
|
|
|
|
func cosh*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the hyperbolic cosine of `z`.
|
|
result = T(0.5) * (exp(z) + exp(-z))
|
|
|
|
func arccosh*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse hyperbolic cosine of `z`.
|
|
result = ln(z + sqrt(z*z - T(1.0)))
|
|
|
|
func tanh*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the hyperbolic tangent of `z`.
|
|
result = sinh(z) / cosh(z)
|
|
|
|
func arctanh*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse hyperbolic tangent of `z`.
|
|
result = T(0.5) * (ln((T(1.0)+z) / (T(1.0)-z)))
|
|
|
|
func coth*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the hyperbolic cotangent of `z`.
|
|
result = cosh(z) / sinh(z)
|
|
|
|
func arccoth*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse hyperbolic cotangent of `z`.
|
|
result = T(0.5) * (ln(T(1.0) + T(1.0)/z) - ln(T(1.0) - T(1.0)/z))
|
|
|
|
func sech*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the hyperbolic secant of `z`.
|
|
result = T(2.0) / (exp(z) + exp(-z))
|
|
|
|
func arcsech*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse hyperbolic secant of `z`.
|
|
result = ln(1.0/z + sqrt(T(1.0)/z+T(1.0)) * sqrt(T(1.0)/z-T(1.0)))
|
|
|
|
func csch*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the hyperbolic cosecant of `z`.
|
|
result = T(2.0) / (exp(z) - exp(-z))
|
|
|
|
func arccsch*[T](z: Complex[T]): Complex[T] =
|
|
## Returns the inverse hyperbolic cosecant of `z`.
|
|
result = ln(T(1.0)/z + sqrt(T(1.0)/(z*z) + T(1.0)))
|
|
|
|
func phase*[T](z: Complex[T]): T =
|
|
## Returns the phase (or argument) of `z`, that is the angle in polar representation.
|
|
##
|
|
## | `result = arctan2(z.im, z.re)`
|
|
arctan2(z.im, z.re)
|
|
|
|
func polar*[T](z: Complex[T]): tuple[r, phi: T] =
|
|
## Returns `z` in polar coordinates.
|
|
##
|
|
## | `result.r = abs(z)`
|
|
## | `result.phi = phase(z)`
|
|
##
|
|
## **See also:**
|
|
## * `rect func<#rect,T,T>`_ for the inverse operation
|
|
(r: abs(z), phi: phase(z))
|
|
|
|
func rect*[T](r, phi: T): Complex[T] =
|
|
## Returns the complex number with polar coordinates `r` and `phi`.
|
|
##
|
|
## | `result.re = r * cos(phi)`
|
|
## | `result.im = r * sin(phi)`
|
|
##
|
|
## **See also:**
|
|
## * `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)"`.
|
|
runnableExamples:
|
|
doAssert $complex(1.0, 2.0) == "(1.0, 2.0)"
|
|
|
|
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.}
|