diff options
Diffstat (limited to 'lib/pure/rationals.nim')
| -rw-r--r-- | lib/pure/rationals.nim | 365 |
1 files changed, 215 insertions, 150 deletions
diff --git a/lib/pure/rationals.nim b/lib/pure/rationals.nim index 04aa8316a..5f806bd70 100644 --- a/lib/pure/rationals.nim +++ b/lib/pure/rationals.nim @@ -8,148 +8,213 @@ # -## This module implements rational numbers, consisting of a numerator `num` and -## a denominator `den`, both of type int. The denominator can not be 0. +## This module implements rational numbers, consisting of a numerator and +## a denominator. The denominator can not be 0. -import math +runnableExamples: + let + r1 = 1 // 2 + r2 = -3 // 4 + + doAssert r1 + r2 == -1 // 4 + doAssert r1 - r2 == 5 // 4 + doAssert r1 * r2 == -3 // 8 + doAssert r1 / r2 == -2 // 3 + +import std/[math, hashes] +when defined(nimPreviewSlimSystem): + import std/assertions type Rational*[T] = object - ## a rational number, consisting of a numerator and denominator + ## A rational number, consisting of a numerator `num` and a denominator `den`. num*, den*: T -proc initRational*[T](num, den: T): Rational[T] = - ## Create a new rational number. +func reduce*[T: SomeInteger](x: var Rational[T]) = + ## Reduces the rational number `x`, so that the numerator and denominator + ## have no common divisors other than 1 (and -1). + ## If `x` is 0, raises `DivByZeroDefect`. + ## + ## **Note:** This is called automatically by the various operations on rationals. + runnableExamples: + var r = Rational[int](num: 2, den: 4) # 1/2 + reduce(r) + doAssert r.num == 1 + doAssert r.den == 2 + if x.den == 0: + raise newException(DivByZeroDefect, "division by zero") + let common = gcd(x.num, x.den) + if x.den > 0: + x.num = x.num div common + x.den = x.den div common + when T isnot SomeUnsignedInt: + if x.den < 0: + x.num = -x.num div common + x.den = -x.den div common + +func initRational*[T: SomeInteger](num, den: T): Rational[T] = + ## Creates a new rational number with numerator `num` and denominator `den`. + ## `den` must not be 0. + ## + ## **Note:** `den != 0` is not checked when assertions are turned off. + assert(den != 0, "a denominator of zero is invalid") result.num = num result.den = den + reduce(result) -proc `//`*[T](num, den: T): Rational[T] = initRational[T](num, den) - ## A friendlier version of `initRational`. Example usage: - ## - ## .. code-block:: nim - ## var x = 1//3 + 1//5 +func `//`*[T](num, den: T): Rational[T] = + ## A friendlier version of `initRational <#initRational,T,T>`_. + runnableExamples: + let x = 1 // 3 + 1 // 5 + doAssert x == 8 // 15 + + initRational[T](num, den) + +func `$`*[T](x: Rational[T]): string = + ## Turns a rational number into a string. + runnableExamples: + doAssert $(1 // 2) == "1/2" -proc `$`*[T](x: Rational[T]): string = - ## Turn a rational number into a string. result = $x.num & "/" & $x.den -proc toRational*[T](x: SomeInteger): Rational[T] = - ## Convert some integer `x` to a rational number. +func toRational*[T: SomeInteger](x: T): Rational[T] = + ## Converts some integer `x` to a rational number. + runnableExamples: + doAssert toRational(42) == 42 // 1 + result.num = x result.den = 1 -proc toFloat*[T](x: Rational[T]): float = - ## Convert a rational number `x` to a float. +func toRational*(x: float, + n: int = high(int) shr (sizeof(int) div 2 * 8)): Rational[int] = + ## Calculates the best rational approximation of `x`, + ## where the denominator is smaller than `n` + ## (default is the largest possible `int` for maximal resolution). + ## + ## The algorithm is based on the theory of continued fractions. + # David Eppstein / UC Irvine / 8 Aug 1993 + # With corrections from Arno Formella, May 2008 + runnableExamples: + let x = 1.2 + doAssert x.toRational.toFloat == x + + var + m11, m22 = 1 + m12, m21 = 0 + ai = int(x) + x = x + while m21 * ai + m22 <= n: + swap m12, m11 + swap m22, m21 + m11 = m12 * ai + m11 + m21 = m22 * ai + m21 + if x == float(ai): break # division by zero + x = 1 / (x - float(ai)) + if x > float(high(int32)): break # representation failure + ai = int(x) + result = m11 // m21 + +func toFloat*[T](x: Rational[T]): float = + ## Converts 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 +func toInt*[T](x: Rational[T]): int = + ## Converts 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. +func `+`*[T](x, y: Rational[T]): Rational[T] = + ## Adds 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`. +func `+`*[T](x: Rational[T], y: T): Rational[T] = + ## Adds the rational `x` to the 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`. +func `+`*[T](x: T, y: Rational[T]): Rational[T] = + ## Adds the int `x` to the 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`. +func `+=`*[T](x: var Rational[T], y: Rational[T]) = + ## Adds the rational `y` to the rational `x` in-place. 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`. +func `+=`*[T](x: var Rational[T], y: T) = + ## Adds the int `y` to the rational `x` in-place. x.num += y * x.den -proc `-` *[T](x: Rational[T]): Rational[T] = +func `-`*[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. +func `-`*[T](x, y: Rational[T]): Rational[T] = + ## Subtracts 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`. +func `-`*[T](x: Rational[T], y: T): Rational[T] = + ## Subtracts the int `y` from the 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 +func `-`*[T](x: T, y: Rational[T]): Rational[T] = + ## Subtracts the rational `y` from the 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`. +func `-=`*[T](x: var Rational[T], y: Rational[T]) = + ## Subtracts the rational `y` from the rational `x` in-place. 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`. +func `-=`*[T](x: var Rational[T], y: T) = + ## Subtracts the int `y` from the rational `x` in-place. x.num -= y * x.den -proc `*` *[T](x, y: Rational[T]): Rational[T] = - ## Multiply two rational numbers. +func `*`*[T](x, y: Rational[T]): Rational[T] = + ## Multiplies 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`. +func `*`*[T](x: Rational[T], y: T): Rational[T] = + ## Multiplies the rational `x` with the 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`. +func `*`*[T](x: T, y: Rational[T]): Rational[T] = + ## Multiplies the int `x` with the 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`. +func `*=`*[T](x: var Rational[T], y: Rational[T]) = + ## Multiplies the rational `x` by `y` in-place. x.num *= y.num x.den *= y.den reduce(x) -proc `*=` *[T](x: var Rational[T], y: T) = - ## Multiply int `y` to rational `x`. +func `*=`*[T](x: var Rational[T], y: T) = + ## Multiplies the rational `x` by the int `y` in-place. x.num *= y reduce(x) -proc reciprocal*[T](x: Rational[T]): Rational[T] = - ## Calculate the reciprocal of `x`. (1/x) +func reciprocal*[T](x: Rational[T]): Rational[T] = + ## Calculates the reciprocal of `x` (`1/x`). + ## If `x` is 0, raises `DivByZeroDefect`. if x.num > 0: result.num = x.den result.den = x.num @@ -157,119 +222,119 @@ proc reciprocal*[T](x: Rational[T]): Rational[T] = result.num = -x.den result.den = -x.num else: - raise newException(DivByZeroError, "division by zero") + raise newException(DivByZeroDefect, "division by zero") -proc `/`*[T](x, y: Rational[T]): Rational[T] = - ## Divide rationals `x` by `y`. +func `/`*[T](x, y: Rational[T]): Rational[T] = + ## Divides the rational `x` by the rational `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`. +func `/`*[T](x: Rational[T], y: T): Rational[T] = + ## Divides the rational `x` by the 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`. +func `/`*[T](x: T, y: Rational[T]): Rational[T] = + ## Divides the int `x` by the 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. +func `/=`*[T](x: var Rational[T], y: Rational[T]) = + ## Divides the rational `x` by the rational `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. +func `/=`*[T](x: var Rational[T], y: T) = + ## Divides the rational `x` by the int `y` in-place. x.den *= y reduce(x) -proc cmp*(x, y: Rational): int = - ## Compares two rationals. +func cmp*(x, y: Rational): int = + ## Compares two rationals. Returns + ## * a value less than zero, if `x < y` + ## * a value greater than zero, if `x > y` + ## * zero, if `x == y` (x - y).num -proc `<` *(x, y: Rational): bool = +func `<`*(x, y: Rational): bool = + ## Returns true if `x` is less than `y`. (x - y).num < 0 -proc `<=` *(x, y: Rational): bool = +func `<=`*(x, y: Rational): bool = + ## Returns tue if `x` is less than or equal to `y`. (x - y).num <= 0 -proc `==` *(x, y: Rational): bool = +func `==`*(x, y: Rational): bool = + ## Compares two rationals for equality. (x - y).num == 0 -proc abs*[T](x: Rational[T]): Rational[T] = +func abs*[T](x: Rational[T]): Rational[T] = + ## Returns the absolute value of `x`. + runnableExamples: + doAssert abs(1 // 2) == 1 // 2 + doAssert abs(-1 // 2) == 1 // 2 + 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 = -1 // -2 - m1 = -1 // 1 - tt = 10 // 2 - - assert( a == a ) - assert( (a-a) == z ) - assert( (a+b) == o ) - assert( (a/b) == o ) - assert( (a*b) == 1 // 4 ) - assert( (3/a) == 6 // 1 ) - assert( (a/3) == 1 // 6 ) - assert( a*b == 1 // 4 ) - assert( tt*z == z ) - assert( 10*a == tt ) - assert( a*10 == tt ) - assert( tt/10 == a ) - assert( a-m1 == 3 // 2 ) - assert( a+m1 == -1 // 2 ) - assert( m1+tt == 16 // 4 ) - assert( m1-tt == 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 = 1//3 - - x *= 5//1 - assert( x == 5//3 ) - x += 2 // 9 - assert( x == 17//9 ) - x -= 9//18 - assert( x == 25//18 ) - x /= 1//2 - assert( x == 50//18 ) - - var y = 1//3 - - y *= 4 - assert( y == 4//3 ) - y += 5 - assert( y == 19//3 ) - y -= 2 - assert( y == 13//3 ) - y /= 9 - assert( y == 13//27 ) - - assert toRational[int, int](5) == 5//1 - assert abs(toFloat(y) - 0.4814814814814815) < 1.0e-7 - assert toInt(z) == 0 +func `div`*[T: SomeInteger](x, y: Rational[T]): T = + ## Computes the rational truncated division. + (x.num * y.den) div (y.num * x.den) + +func `mod`*[T: SomeInteger](x, y: Rational[T]): Rational[T] = + ## Computes the rational modulo by truncated division (remainder). + ## This is same as `x - (x div y) * y`. + result = ((x.num * y.den) mod (y.num * x.den)) // (x.den * y.den) + reduce(result) + +func floorDiv*[T: SomeInteger](x, y: Rational[T]): T = + ## Computes the rational floor division. + ## + ## Floor division is conceptually defined as `floor(x / y)`. + ## This is different from the `div` operator, which is defined + ## as `trunc(x / y)`. That is, `div` rounds towards 0 and `floorDiv` + ## rounds down. + floorDiv(x.num * y.den, y.num * x.den) + +func floorMod*[T: SomeInteger](x, y: Rational[T]): Rational[T] = + ## Computes the rational modulo by floor division (modulo). + ## + ## This is same as `x - floorDiv(x, y) * y`. + ## This func behaves the same as the `%` operator in Python. + result = floorMod(x.num * y.den, y.num * x.den) // (x.den * y.den) + reduce(result) + +func hash*[T](x: Rational[T]): Hash = + ## Computes the hash for the rational `x`. + # reduce first so that hash(x) == hash(y) for x == y + var copy = x + reduce(copy) + + var h: Hash = 0 + h = h !& hash(copy.num) + h = h !& hash(copy.den) + result = !$h + +func `^`*[T: SomeInteger](x: Rational[T], y: T): Rational[T] = + ## Computes `x` to the power of `y`. + ## + ## The exponent `y` must be an integer. Negative exponents are supported + ## but floating point exponents are not. + runnableExamples: + doAssert (-3 // 5) ^ 0 == (1 // 1) + doAssert (-3 // 5) ^ 1 == (-3 // 5) + doAssert (-3 // 5) ^ 2 == (9 // 25) + doAssert (-3 // 5) ^ -2 == (25 // 9) + + if y >= 0: + result.num = x.num ^ y + result.den = x.den ^ y + else: + result.num = x.den ^ -y + result.den = x.num ^ -y + # Note that all powers of reduced rationals are already reduced, + # so we don't need to call reduce() here |