Merge pull request #1840 from def-/extend-math

Rational numbers and a few additions to math and complex

doc/lib.txt

@@ -185,6 +185,9 @@ Math libraries
 * `complex <complex.html>`_
   This module implements complex numbers and their mathematical operations.
 
+* `rationals <rationals.html>`_
+  This module implements rational numbers and their mathematical operations.
+
 * `fenv <fenv.html>`_
   Floating-point environment. Handling of floating-point rounding and
   exceptions (overflow, zero-devide, etc.).

lib/pure/complex.nim

@@ -27,6 +27,11 @@ type
 
 {.deprecated: [TComplex: Complex].}
 
+proc toComplex*(x: SomeInteger): Complex =
+  ## Convert some integer ``x`` to a complex number.
+  result.re = x
+  result.im = 0
+
 proc `==` *(x, y: Complex): bool =
   ## Compare two complex numbers `x` and `y` for equality.
   result = x.re == y.re and x.im == y.im
@@ -291,6 +296,21 @@ proc cosh*(z: Complex): Complex =
   result = 0.5*(exp(z)+exp(-z))
 
 
+proc phase*(z: Complex): float =
+  ## Returns the phase of `z`.
+  arctan2(z.im, z.re)
+
+proc polar*(z: Complex): tuple[r, phi: float] =
+  ## Returns `z` in polar coordinates.
+  result.r = abs(z)
+  result.phi = phase(z)
+
+proc rect*(r: float, phi: float): Complex =
+  ## Returns the complex number with poolar coordinates `r` and `phi`.
+  result.re = r * cos(phi)
+  result.im = sin(phi)
+
+
 proc `$`*(z: Complex): string =
   ## Returns `z`'s string representation as ``"(re, im)"``.
   result = "(" & $z.re & ", " & $z.im & ")"
@@ -344,6 +364,9 @@ when isMainModule:
   assert( arcsin(a) =~ (0.427078586392476, 1.528570919480998) )
   assert( arccos(a) =~ (1.14371774040242, -1.52857091948100) )
 
-  assert( cosh(a) =~ (-0.642148124715520, 1.068607421382778) ) 
+  assert( cosh(a) =~ (-0.642148124715520, 1.068607421382778) )
   assert( sinh(a) =~ (-0.489056259041294, 1.403119250622040) )
-  
+
+  assert( phase(a) == 1.1071487177940904 )
+  assert( polar(a) =~ (2.23606797749979, 1.1071487177940904) )
+  assert( rect(1.0, 2.0) =~ (-0.4161468365471424, 0.9092974268256817) )

lib/pure/math.nim

@@ -280,7 +280,7 @@ proc random*[T](x: Slice[T]): T =
   ## For a slice `a .. b` returns a value in the range `a .. b-1`.
   result = random(x.b - x.a) + x.a
 
-proc random[T](a: openArray[T]): T =
+proc random*[T](a: openArray[T]): T =
   ## returns a random element from the openarray `a`.
   result = a[random(a.low..a.len)]
 
@@ -329,6 +329,31 @@ proc standardDeviation*(s: RunningStat): float =
 {.pop.}
 {.pop.}
 
+proc `^`*[T](x, y: T): T =
+  ## Computes ``x`` to the power ``y`. ``x`` must be non-negative, use
+  ## `pow <#pow,float,float>` for negative exponents.
+  assert y >= 0
+  var (x, y) = (x, y)
+  result = 1
+
+  while y != 0:
+    if (y and 1) != 0:
+      result *= x
+    y = y shr 1
+    x *= x
+
+proc gcd*[T](x, y: T): T =
+  ## Computes the greatest common divisor of ``x`` and ``y``.
+  var (x,y) = (x,y)
+  while y != 0:
+    x = x mod y
+    swap x, y
+  abs x
+
+proc lcm*[T](x, y: T): T =
+  ## Computes the least common multiple of ``x`` and ``y``.
+  x div gcd(x, y) * y
+
 when isMainModule and not defined(JS):
   proc gettime(dummy: ptr cint): cint {.importc: "time", header: "<time.h>".}
 

lib/pure/rationals.nim

@@ -0,0 +1,265 @@
+#
+#
+#            Nim's Runtime Library
+#        (c) Copyright 2015 Dennis Felsing
+#
+#    See the file "copying.txt", included in this
+#    distribution, for details about the copyright.
+#
+
+
+## This module implements rational numbers, consisting of a numerator `num` and
+## a denominator `den`, both of type int. The denominator can not be 0.
+
+import math
+
+type Rational*[T] = object
+  ## a rational number, consisting of a numerator and denominator
+  num*, den*: T
+
+proc initRational*[T](num, den: T): Rational[T] =
+  ## Create a new rational number.
+  result.num = num
+  result.den = den
+
+proc toRational*[T](x: SomeInteger): Rational[T] =
+  ## Convert some integer `x` to a rational number.
+  result.num = x
+  result.den = 1
+
+proc toFloat*[T](x: Rational[T]): float =
+  ## Convert a rational number `x` to a float.
+  x.num / x.den
+
+proc toInt*[T](x: Rational[T]): int =
+  ## Convert a rational number `x` to an int. Conversion rounds towards 0 if
+  ## `x` does not contain an integer value.
+  x.num div x.den
+
+proc reduce*[T](x: var Rational[T]) =
+  ## Reduce rational `x`.
+  let common = gcd(x.num, x.den)
+  if x.den > 0:
+    x.num = x.num div common
+    x.den = x.den div common
+  elif x.den < 0:
+    x.num = -x.num div common
+    x.den = -x.den div common
+  else:
+    raise newException(DivByZeroError, "division by zero")
+
+proc `+` *[T](x, y: Rational[T]): Rational[T] =
+  ## Add two rational numbers.
+  let common = lcm(x.den, y.den)
+  result.num = common div x.den * x.num + common div y.den * y.num
+  result.den = common
+  reduce(result)
+
+proc `+` *[T](x: Rational[T], y: T): Rational[T] =
+  ## Add rational `x` to int `y`.
+  result.num = x.num + y * x.den
+  result.den = x.den
+
+proc `+` *[T](x: T, y: Rational[T]): Rational[T] =
+  ## Add int `x` to rational `y`.
+  result.num = x * y.den + y.num
+  result.den = y.den
+
+proc `+=` *[T](x: var Rational[T], y: Rational[T]) =
+  ## Add rational `y` to rational `x`.
+  let common = lcm(x.den, y.den)
+  x.num = common div x.den * x.num + common div y.den * y.num
+  x.den = common
+  reduce(x)
+
+proc `+=` *[T](x: var Rational[T], y: T) =
+  ## Add int `y` to rational `x`.
+  x.num += y * x.den
+
+proc `-` *[T](x: Rational[T]): Rational[T] =
+  ## Unary minus for rational numbers.
+  result.num = -x.num
+  result.den = x.den
+
+proc `-` *[T](x, y: Rational[T]): Rational[T] =
+  ## Subtract two rational numbers.
+  let common = lcm(x.den, y.den)
+  result.num = common div x.den * x.num - common div y.den * y.num
+  result.den = common
+  reduce(result)
+
+proc `-` *[T](x: Rational[T], y: T): Rational[T] =
+  ## Subtract int `y` from rational `x`.
+  result.num = x.num - y * x.den
+  result.den = x.den
+
+proc `-` *[T](x: T, y: Rational[T]): Rational[T] =
+  ## Subtract rational `y` from int `x`.
+  result.num = - x * y.den + y.num
+  result.den = y.den
+
+proc `-=` *[T](x: var Rational[T], y: Rational[T]) =
+  ## Subtract rational `y` from rational `x`.
+  let common = lcm(x.den, y.den)
+  x.num = common div x.den * x.num - common div y.den * y.num
+  x.den = common
+  reduce(x)
+
+proc `-=` *[T](x: var Rational[T], y: T) =
+  ## Subtract int `y` from rational `x`.
+  x.num -= y * x.den
+
+proc `*` *[T](x, y: Rational[T]): Rational[T] =
+  ## Multiply two rational numbers.
+  result.num = x.num * y.num
+  result.den = x.den * y.den
+  reduce(result)
+
+proc `*` *[T](x: Rational[T], y: T): Rational[T] =
+  ## Multiply rational `x` with int `y`.
+  result.num = x.num * y
+  result.den = x.den
+  reduce(result)
+
+proc `*` *[T](x: T, y: Rational[T]): Rational[T] =
+  ## Multiply int `x` with rational `y`.
+  result.num = x * y.num
+  result.den = y.den
+  reduce(result)
+
+proc `*=` *[T](x: var Rational[T], y: Rational[T]) =
+  ## Multiply rationals `y` to `x`.
+  x.num *= y.num
+  x.den *= y.den
+  reduce(x)
+
+proc `*=` *[T](x: var Rational[T], y: T) =
+  ## Multiply int `y` to rational `x`.
+  x.num *= y
+  reduce(x)
+
+proc reciprocal*[T](x: Rational[T]): Rational[T] =
+  ## Calculate the reciprocal of `x`. (1/x)
+  if x.num > 0:
+    result.num = x.den
+    result.den = x.num
+  elif x.num < 0:
+    result.num = -x.den
+    result.den = -x.num
+  else:
+    raise newException(DivByZeroError, "division by zero")
+
+proc `/`*[T](x, y: Rational[T]): Rational[T] =
+  ## Divide rationals `x` by `y`.
+  result.num = x.num * y.den
+  result.den = x.den * y.num
+  reduce(result)
+
+proc `/`*[T](x: Rational[T], y: T): Rational[T] =
+  ## Divide rational `x` by int `y`.
+  result.num = x.num
+  result.den = x.den * y
+  reduce(result)
+
+proc `/`*[T](x: T, y: Rational[T]): Rational[T] =
+  ## Divide int `x` by Rational `y`.
+  result.num = x * y.den
+  result.den = y.num
+  reduce(result)
+
+proc `/=`*[T](x: var Rational[T], y: Rational[T]) =
+  ## Divide rationals `x` by `y` in place.
+  x.num *= y.den
+  x.den *= y.num
+  reduce(x)
+
+proc `/=`*[T](x: var Rational[T], y: T) =
+  ## Divide rational `x` by int `y` in place.
+  x.den *= y
+  reduce(x)
+
+proc cmp*(x, y: Rational): int =
+  ## Compares two rationals.
+  (x - y).num
+
+proc `<` *(x, y: Rational): bool =
+  (x - y).num < 0
+
+proc `<=` *(x, y: Rational): bool =
+  (x - y).num <= 0
+
+proc `==` *(x, y: Rational): bool =
+  (x - y).num == 0
+
+proc abs*[T](x: Rational[T]): Rational[T] =
+  result.num = abs x.num
+  result.den = abs x.den
+
+when isMainModule:
+  var
+    z = Rational[int](num: 0, den: 1)
+    o = initRational(num=1, den=1)
+    a = initRational(1, 2)
+    b = initRational(-1, -2)
+    m1 = initRational(-1, 1)
+    tt = initRational(10, 2)
+
+  assert( a     == a )
+  assert( (a-a) == z )
+  assert( (a+b) == o )
+  assert( (a/b) == o )
+  assert( (a*b) == initRational(1, 4) )
+  assert( (3/a) == initRational(6,1) )
+  assert( (a/3) == initRational(1,6) )
+  assert( a*b   == initRational(1,4) )
+  assert( tt*z  == z )
+  assert( 10*a  == tt )
+  assert( a*10  == tt )
+  assert( tt/10 == a  )
+  assert( a-m1  == initRational(3, 2) )
+  assert( a+m1  == initRational(-1, 2) )
+  assert( m1+tt == initRational(16, 4) )
+  assert( m1-tt == initRational(6, -1) )
+
+  assert( z < o )
+  assert( z <= o )
+  assert( z == z )
+  assert( cmp(z, o) < 0 )
+  assert( cmp(o, z) > 0 )
+
+  assert( o == o )
+  assert( o >= o )
+  assert( not(o > o) )
+  assert( cmp(o, o) == 0 )
+  assert( cmp(z, z) == 0 )
+
+  assert( a == b )
+  assert( a >= b )
+  assert( not(b > a) )
+  assert( cmp(a, b) == 0 )
+
+  var x = initRational(1,3)
+
+  x *= initRational(5,1)
+  assert( x == initRational(5,3) )
+  x += initRational(2,9)
+  assert( x == initRational(17,9) )
+  x -= initRational(9,18)
+  assert( x == initRational(25,18) )
+  x /= initRational(1,2)
+  assert( x == initRational(50,18) )
+
+  var y = initRational(1,3)
+
+  y *= 4
+  assert( y == initRational(4,3) )
+  y += 5
+  assert( y == initRational(19,3) )
+  y -= 2
+  assert( y == initRational(13,3) )
+  y /= 9
+  assert( y == initRational(13,27) )
+
+  assert toRational[int, int](5) == initRational(5,1)
+  assert abs(toFloat(y) - 0.4814814814814815) < 1.0e-7
+  assert toInt(z) == 0