diff options
| -rw-r--r-- | doc/lib.txt | 3 | ||||
| -rw-r--r-- | lib/pure/complex.nim | 27 | ||||
| -rw-r--r-- | lib/pure/math.nim | 27 | ||||
| -rw-r--r-- | lib/pure/rationals.nim | 265 |
4 files changed, 319 insertions, 3 deletions
diff --git a/doc/lib.txt b/doc/lib.txt index 1f972179d..76920c6a9 100644 --- a/doc/lib.txt +++ b/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.). diff --git a/lib/pure/complex.nim b/lib/pure/complex.nim index 04cffe4e4..9f1546eed 100644 --- a/lib/pure/complex.nim +++ b/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) ) - \ No newline at end of file + + assert( phase(a) == 1.1071487177940904 ) + assert( polar(a) =~ (2.23606797749979, 1.1071487177940904) ) + assert( rect(1.0, 2.0) =~ (-0.4161468365471424, 0.9092974268256817) ) diff --git a/lib/pure/math.nim b/lib/pure/math.nim index b25a1df3a..c902af381 100644 --- a/lib/pure/math.nim +++ b/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>".} diff --git a/lib/pure/rationals.nim b/lib/pure/rationals.nim new file mode 100644 index 000000000..40c61f1d9 --- /dev/null +++ b/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 |