added rational numbers

src/mmath/hmatrix.nim

@@ -42,7 +42,7 @@ proc newHMatrix*[T](rows, columns : int, values : openarray[T]) : HMatrix[T] =
 proc identity*[T](sz : int) : HMatrix[T] =
     result = newHMatrix[T](sz,sz,0)
     for i in 0...sz:
-        result[i,i] = 1
+        result[i,i] = T.one
 
 proc `$`*[T](m : HMatrix[T]) : string = 
     result = "["
@@ -60,7 +60,7 @@ proc `$`*[T](m : HMatrix[T]) : string =
         else: result &= ", "
 
 proc `*`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] =
-    result = newHMatrix[T](m1.rows, m2.columns, 0)
+    result = newHMatrix[T](m1.rows, m2.columns, T.zero)
     for i in 0...result.rows:
         for j in 0...result.columns:
             for k in 0...m1.columns:
@@ -239,7 +239,7 @@ proc det*[T](m : HMatrix[T]) : T =
         raise newException(SizeError, "Matrix must be square in order to compute the determinant")
     var clone = m.clone()
     clone.gauss_jordan_low()
-    result = 1
+    result = T.one
     for i in 0...clone.rows:
         result *= clone[i, i]
 
@@ -256,7 +256,7 @@ proc invert*[T](m : HMatrix[T]) : HMatrix[T] =
             result[i, j] /= f
 
 proc triu*[T](m : HMatrix[T], diag_replace: T) : HMatrix[T] =
-    result = newHMatrix[T](m.rows, m.columns, 0)
+    result = newHMatrix[T](m.rows, m.columns, T.zero)
     for i in 0...m.rows:
         for j in i...m.columns:
             if i == j:
@@ -265,14 +265,14 @@ proc triu*[T](m : HMatrix[T], diag_replace: T) : HMatrix[T] =
                 result[i, j] = m[i, j]
 
 proc triu*[T](m : HMatrix[T]) : HMatrix[T] =
-    result = newHMatrix[T](m.rows, m.columns, 0)
+    result = newHMatrix[T](m.rows, m.columns, T.zero)
     for i in 0...m.rows:
         for j in i...m.columns:
             result[i, j] = m[i, j]
                 
                 
 proc tril*[T](m : HMatrix[T], diag_replacement : T) : HMatrix[T] =
-    result = newHMatrix[T](m.rows, m.columns, 0)
+    result = newHMatrix[T](m.rows, m.columns, T.zero)
     for i in 0...m.rows:
         for j in 0...(i + 1):
             if i == j:
@@ -281,13 +281,13 @@ proc tril*[T](m : HMatrix[T], diag_replacement : T) : HMatrix[T] =
                 result[i, j] = m[i, j]
 
 proc tril*[T](m : HMatrix[T]) : HMatrix[T] =
-    result = newHMatrix[T](m.rows, m.columns, 0)
+    result = newHMatrix[T](m.rows, m.columns, T.zero)
     for i in 0...m.rows:
         for j in 0...(i + 1):
             result[i, j] = m[i, j]
 
 proc lu_row[T](m : var HMatrix[T], i : int) =
-    if m[i, i] == 0:
+    if m[i, i] == T.zero:
         raise newException(SingularMatrixError, "Matrix is singular")
     for j in i...m.columns:
         for k in 0...i:
@@ -345,9 +345,9 @@ proc lu_invert*[T](m : HMatrix[T], pivot : HPivot[T]) : HMatrix[T] =
     for i in 0...m.rows:
         for j in 0...m.rows:
             if pivot[j] == i:
-                result[j, i] = 1
+                result[j, i] = T.one
             else:
-                result[j, i] = 0
+                result[j, i] = T.zero
             for k in range(j):
                 result[j, i] -= m[j, k] * result[k, i]
         for j in m.rows-->0:
@@ -361,11 +361,11 @@ proc lu_det*[T](m : var HMatrix[T]) : T =
     if m.rows != m.columns:
         raise newException(SizeError, "Matrix must be square in order to compute the determinant")
     let pivot = m.lup()
-    result = 1
+    result = T.one
     for i in 0...m.rows:
         result *= m[i, i]
     if pivot.permutations mod 2 != 0:
-        result *= -1
+        result *= -T.one
 
 proc lu_det*[T](m : HMatrix[T]) : T = 
     if m.rows != m.columns:
@@ -374,7 +374,7 @@ proc lu_det*[T](m : HMatrix[T]) : T =
     lu_det(clone)
 
 proc squared_norm2*[T](m : HMatrix[T]): T =
-    result = T(0)
+    result = T.zero
     for i, j, v in items(m):
         result += v * v
 
@@ -394,5 +394,5 @@ proc `*`*[T](pivot : HPivot[T], m : HMatrix[T]) : HMatrix[T] =
 proc from_pivot*[T](pivot : HPivot[T]): HMatrix[T] =
     result = newHMatrix[T](len(pivot), len(pivot))
     for i in 0...pivot.len:
-        result[pivot[i], i] = 1
+        result[pivot[i], i] = T.one
 

src/mmath/rational.nim

@@ -0,0 +1,118 @@
+proc one*[T : SomeNumber](_: typedesc[T]) : T = 1.T
+
+proc zero*[T : SomeNumber](_: typedesc[T]) : T = 0.T
+
+proc gcd*[T : SomeInteger](a : T, b : T) : T =
+    var n1 = a
+    var n2 = b 
+    var tmp : T
+    while n2 != T.zero:
+        tmp = n1
+        n1 = n2
+        n2 = tmp mod n2
+    return n1
+
+proc mcm*[T : SomeInteger](n1 : T, n2 : T) : T = n1 * n2 div gcd(n1, n2)
+
+proc pow*[T : SomeInteger](n : T, e : T) : T =
+    result = T.one
+    for i in 0..<e:
+        result *= n    
+
+type Rational*[T : SomeInteger] = object
+    num* : T
+    den* : T
+
+proc simplify*[T](self : var Rational[T]) : void =
+    let gcd = gcd(self.num.abs(), self.den.abs())
+    let invertSign = self.den < T.zero
+    let num = self.num div gcd * (if invertSign: (T.zero - T.one) else : T.one)
+    let den = self.den div gcd * (if invertSign: (T.zero - T.one) else : T.one)
+    self.num = num
+    self.den = den
+
+proc newRational*[T](num : T) : Rational[T] = Rational[T](num: num, den: T(1))
+proc newRational*[T](num : T, den: T) : Rational[T] = Rational[T](num: num, den: den)
+
+proc zero*[T](_: typedesc[Rational[T]]): Rational[T] = newRational(T.zero, T.one)
+
+proc one*[T](_: typedesc[Rational[T]]): Rational[T] = newRational(T.one, T.one)
+
+proc abs*[T](self : Rational[T]): Rational[T] = Rational[T](num: self.num.abs(), den: self.den.abs())
+
+proc pow*[T](self : Rational[T], e : SomeInteger): Rational[T] = newRational(self.num.pow(e), self.den.pow(e))
+
+proc sqrt*[T](self : Rational[T]) : Rational[T] = newRational(sqrt(self.num), sqrt(self.den))
+
+proc `-`*[T](r1 : Rational[T]) : Rational[T] = newRational(-r1.num, r1.den)
+
+proc `+`*[T](r1 : Rational[T], r2: Rational[T]) : Rational[T] =
+    let den = mcm(r1.den, r2.den)
+    result = newRational(r1.num * den div r1.den + r2.num * den div r2.den, den)
+    result.simplify()
+
+proc `-`*[T](r1 : Rational[T], r2: Rational[T]) : Rational[T] =
+    let den = mcm(r1.den, r2.den)
+    result = newRational(r1.num * den div r1.den - r2.num * den div r2.den, den)
+    result.simplify()
+
+proc `*`*[T](r1 : Rational[T], r2: Rational[T]) : Rational[T] = 
+    result = newRational(r1.num * r2.num, r1.den * r2.den)
+    result.simplify()
+
+proc `/`*[T](r1 : Rational[T], r2: Rational[T]) : Rational[T] = 
+    result = newRational(r1.num * r2.den, r1.den * r2.num)
+    result.simplify()
+
+proc `+=`*[T](r1 : var Rational[T], r2: Rational[T]) : void =
+    let den = mcm(r1.den, r2.den)
+    r1.num = r1.num * den div r1.den + r2.num * den div r2.den 
+    r1.den = den
+    r1.simplify()
+
+proc `-=`*[T](r1 : var Rational[T], r2: Rational[T]) : void =
+    let den = mcm(r1.den, r2.den)
+    r1.num = r1.num * den div r1.den - r2.num * den div r2.den 
+    r1.den = den
+    r1.simplify()
+
+proc `*=`*[T](r1 : var Rational[T], r2: Rational[T]) : void =
+    r1.num = r1.num * r2.num
+    r1.den = r1.den * r2.den
+    r1.simplify()
+
+proc `/=`*[T](r1 : var Rational[T], r2: Rational[T]) : void =
+    r1.num = r1.num * r2.den 
+    r1.den = r1.den * r2.num
+    r1.simplify()   
+
+proc `cmp`*[T](r1 : Rational[T], r2 : Rational[T]) : int =
+    cmp(r1.num * r2.den, r1.den * r2.num)
+
+proc `==`*[T](r1 : Rational[T], r2 : Rational[T]) : bool =
+    cmp(r1, r2) == 0
+
+proc `<`*[T](r1 : Rational[T], r2 : Rational[T]) : bool =
+    cmp(r1, r2) < 0
+
+proc `<=`*[T](r1 : Rational[T], r2 : Rational[T]) : bool =
+    cmp(r1, r2) <= 0
+
+proc `>`*[T](r1 : Rational[T], r2 : Rational[T]) : bool =
+    cmp(r1, r2) > 0
+
+proc `>=`*[T](r1 : Rational[T], r2 : Rational[T]) : bool =
+    cmp(r1, r2) >= 0
+
+proc `$`*[T](r : Rational[T]): string = 
+    let zero = T.zero
+    let one = T.one
+    if r.num == T.zero and r.den != zero:
+        result = $zero
+    else:
+        let negative = r.num == r.num.abs() xor r.den == r.den.abs()
+        if r.den.abs() == one:
+            result = $(if negative: (zero - one) * r.num else: r.num)  
+        else:
+            result = (if negative: "-" else: "") & $(r.num.abs()) & "/" & $(r.den.abs())
+

src/mmath/smatrix.nim

@@ -38,7 +38,7 @@ proc newSMatrixFromArray*[ROWS, COLUMNS : static[int], T](values : array[0..(ROW
 
 proc identity*[SIZE: static[int], T]() : SquareSMatrix[SIZE, T] =
     for i in 0...SIZE:
-        result[i,i] = 1
+        result[i,i] = T.one
 
 proc `$`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : string = 
     result = "["
@@ -60,8 +60,9 @@ proc `*`*[ROWS1, COLUMNS2, COMMON : static[int], T](
         m2 : SMatrix[COMMON, COLUMNS2, T]) : SMatrix[ROWS1, COLUMNS2, T] =
     for i in 0...result.rows:
         for j in 0...result.columns:
+            result[i, j] = T.zero
             for k in 0...m1.columns:
-                result[i, j] = result[i, j] + m1[i, k] * m2[k, j]
+                result[i, j] += m1[i, k] * m2[k, j]
 
 proc `*`*[SIZE : static[int], T](v : SVector[SIZE, T], m : SquareSMatrix[SIZE, T]) : SVector[SIZE, T] =
     for i in 0...m.rows:
@@ -125,7 +126,7 @@ proc `==`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], m2 :
     return true
 
 proc squared_norm2*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]): T =
-    result = T(0)
+    result = T.zero
     for i, j, v in items(m):
         result += v * v
 
@@ -182,7 +183,7 @@ proc gauss_jordan_low*[ROWS, COLUMNS : static[int], T](
         other : var SMatrix[ROWS, COLUMNS, T]) =
     var pivot = newSPivot[ROWS, T]()
     for i in 0...m.rows:
-        if m[i, i] == 0:
+        if m[i, i] == T.zero:
             for j in (i + 1)...m.columns:
                 if m[j, i] != 0:
                     m.swap_rows(i, j, pivot, other)
@@ -238,7 +239,7 @@ proc gauss_jordan_high*[ROWS, COLUMNS : static[int], T](
 proc det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T =
     var clone = m.clone()
     clone.gauss_jordan_low()
-    result = 1
+    result = T.one
     for i in 0...clone.rows:
         result *= clone[i, i]
 
@@ -252,34 +253,44 @@ proc invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[
         for j in 0...result.columns:
             result[i, j] /= f
 
-proc triu*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replace: T) : SquareSMatrix[SIZE, T] =
+proc triu*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replacement: T) : SquareSMatrix[SIZE, T] =
     for i in 0...m.rows:
-        for j in i...m.columns:
+        for j in 0...m.columns:
-            if i == j:
+            if i < j:
-                result[i, j] = diag_replace
-            else:
                 result[i, j] = m[i, j]
+            elif i == j:
+                result[i, j] = diag_replacement
+            else:
+                result[i, j] = T.zero
 
 proc triu*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] =
     for i in 0...m.rows:
-        for j in i...m.columns:
+        for j in 0...m.columns:
-            result[i, j] = m[i, j]
+            if i <= j:
+                result[i, j] = m[i, j]
+            else:
+                result[i, j] = T.zero
                 
 proc tril*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replacement : T) : SquareSMatrix[SIZE, T] =
     for i in 0...m.rows:
-        for j in 0...(i + 1):
+        for j in 0...m.columns:
-            if i == j:
+            if i > j:
+                result[i, j] = m[i, j]
+            elif i == j:
                 result[i, j] = diag_replacement
             else:
-                result[i, j] = m[i, j]
+                result[i, j] = T.zero
 
 proc tril*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] =
     for i in 0...m.rows:
-        for j in 0...(i + 1):
+        for j in 0...m.columns:
-            result[i, j] = m[i, j]
+            if i >= j:
+                result[i, j] = m[i, j]
+            else:
+                result[i, j] = T.zero
 
 proc lu_row[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T], i : int) =
-    if m[i, i] == 0:
+    if m[i, i] == T.zero:
         raise newException(SingularMatrixError, "Matrix is singular")
     for j in i...m.columns:
         for k in 0...i:
@@ -334,9 +345,9 @@ proc lu_invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], pivot : SPivo
     for i in 0...m.rows:
         for j in 0...m.rows:
             if pivot[j] == i:
-                result[j, i] = 1
+                result[j, i] = T.one
             else:
-                result[j, i] = 0
+                result[j, i] = T.zero
             for k in range(j):
                 result[j, i] -= m[j, k] * result[k, i]
         for j in m.rows-->0:
@@ -348,11 +359,11 @@ proc lu_invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatr
 
 proc lu_det*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) : T =
     let pivot = m.lup()
-    result = 1
+    result = T.one
     for i in 0...m.rows:
         result *= m[i, i]
     if pivot.permutations mod 2 != 0:
-        result *= -1
+        result *= -T.one
 
 proc lu_det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T = 
     var clone = m.clone()
@@ -371,5 +382,5 @@ proc `*`*[ROWS, COLUMNS : static[int], T](pivot : SPivot[ROWS, T], m : SMatrix[R
 proc from_pivot*[SIZE : static[int], T](pivot : SPivot[SIZE, T]): SquareSMatrix[SIZE, T] =
     result = newSMatrix[T](len(pivot), len(pivot))
     for i in 0...SIZE:
-        result[pivot[i], i] = 1
+        result[pivot[i], i] = T.one
 

tests/test_hmatrix.nim

@@ -1,7 +1,8 @@
 from random import initRand, rand
 from nwo/utils import `...`
+from mmath/rational import newRational, Rational, one, zero, abs, `>`, `<`, `==`, `*`, `+`, `-`, `/`, `/=`, `+=`, `sqrt`
 from mmath/hmatrix import det, lu, from_pivot, newHMatrix, HMatrix, lu_det, invert, identity,
-    clone, tril, triu, lu_solve, `*`, `-`, `+`, `+=`, `-=`, `==`, norm2, transpose, lup
+    clone, tril, triu, lu_solve, `*`, `-`, `+`, `+=`, `-=`, `==`, norm2, transpose, lup, squared_norm2
 from mmath/hvector import buildHVector, newHVector, `-`, abs, norm
 import unittest
 
@@ -58,16 +59,16 @@ suite "Nim linear algebra library":
 
     test "LU decomposition":
         var rng = initRand(101325)
-        var arr : array[0..(25 - 1), float32]
+        var arr : array[0..(25 - 1), Rational[int64]]
         for i in 0..<arr.len:
-            arr[i] = rng.rand(-4f32..4f32)    
+            arr[i] = newRational(rng.rand(-20..20).int64, 20.int64)
-        let mtx = newHMatrix[float32](5, 5, arr)
+        let mtx = newHMatrix[Rational[int64]](5, 5, arr)
         var lu = mtx.clone()        
         let pivot = lu.lup()
-        let l = lu.tril(1.0)
+        let l = lu.tril(Rational[int64].one)
         let u = lu.triu()
         let err = pivot * mtx - (l * u)
-        check(err.norm2() < 1e-5)
+        check(err.squared_norm2() == Rational[int64].zero)
 
     test "Linear system solve":
         var rng = initRand(101325)

tests/test_rational.nim

@@ -0,0 +1,96 @@
+import unittest
+from mmath/rational import newRational, `+`, `-`, `*`, `/`, `$`, simplify, gcd, pow, mcm, `==`, `cmp`, `<`, `>`, `<=`, `>=`, abs
+
+suite "Nim rational numbers":
+
+    test "gcd":
+        let a = 14
+        let b = 21
+        check 7 == gcd(a, b)
+
+    test "mcm":
+        let a = 14
+        let b = 21
+        check 42 == mcm(a, b)
+
+    test "pow":
+        check 625 == pow(5, 4)
+
+    test "==":
+        let r1 = newRational(1, 3)
+        let r2 = newRational(5, 15)        
+        let r3 = newRational(35, 105)
+        check r1 == r2
+        check r1 == r3
+        check r2 == r3
+
+    test "+":
+        let r1 = newRational(1, 3)
+        let r2 = newRational(2, 3)
+        let r3 = r1 + r2
+        check r3 == newRational(1)
+
+    test "-":
+        let r1 = newRational(1, 3)
+        let r2 = newRational(2, 3)
+        let r3 = r1 - r2
+        check r3 == newRational(-1, 3)
+
+    test "*":
+        let r1 = newRational(1, 3)
+        let r2 = newRational(2, 3)
+        let r3 = r1 * r2
+        check r3 == newRational(2, 9)
+
+    test "/":
+        let r1 = newRational(1, 3)
+        let r2 = newRational(2, 3)
+        let r3 = r1 / r2
+        check r3 == newRational(1, 2)
+
+    test "simplify":
+        var r1 = newRational(65, 169)
+        check r1 == newRational(5, 13)
+        r1.simplify()
+        check r1.num == 5
+        check r1.den == 13
+
+    test "pow":
+        let r1 = newRational(65, 169)
+        let r2 = r1.pow(2)
+        check r2 == newRational(25, 169)
+
+    test "cmp":
+        let r1 = newRational(65, 169)
+        let r2 = newRational(5, 13)
+        let r3 = newRational(6, 13)
+        check r1 < r3
+        check r3 > r1
+        check r1 <= r3
+        check r1 <= r2
+        check r3 >= r1
+        check r2 >= r1
+        check r2 == r1
+
+    test "abs":
+        let r1 = newRational(-65, -169)
+        let r2 = newRational(65, -169)
+        check r1 == r1.abs()
+        check r1 == r2.abs()
+        check r1 != r2
+
+    test "i64":
+        let p1 = 32452867i64
+        let p2 = 49979687i64
+        let p3 = 15485867i64
+        let r1 = newRational(p3, p1 * p2)
+        let r2 = newRational(p2, p1 * p3)
+        var r3 = newRational(p2, p2 * p3)
+        check newRational(1i64, p3) == r3
+        r3.simplify()
+        check r3.num == 1i64
+        check r3.den == p3
+        expect(OverflowDefect):
+            discard r1 / r2
+        expect(OverflowDefect):
+            discard r1 + r2

tests/test_smatrix.nim

@@ -2,9 +2,10 @@ from random import initRand, rand
 from nwo/utils import `...`
 from mmath/smatrix import det, lu, lup, lu_det, lup, det, from_pivot, invert, identity, 
     newSMatrix, SMatrix, clone, tril, triu, lu_solve, `$`, `*`, `-`, `+`, `-=`, `+=`, `==`, `*`, norm2, newSMatrixFromArray,
-    gauss_jordan_high, gauss_jordan_low, transpose
+    gauss_jordan_high, gauss_jordan_low, transpose, squared_norm2
 from mmath/svector import buildSVector, `-`, abs, norm, Svector
 from mmath/pivot import SingularMatrixError
+from mmath/rational import newRational, Rational, one, zero, abs, `>`, `<`, `==`, `*`, `+`, `-`, `/`, `/=`, `+=`, `sqrt`
 import unittest
 
 suite "Nim linear algebra library":
@@ -58,19 +59,19 @@ suite "Nim linear algebra library":
     test "Inverse":
         let err = test_matrix * test_matrix.invert() - identity[3, float32]()
         check(err.norm2() < 1e-5)
-
+    
     test "LU decomposition":
         var rng = initRand(101325)
-        var arr : array[0..(25 - 1), float32]
+        var arr : array[0..(25 - 1), Rational[int64]]
         for i in 0..<arr.len:
-            arr[i] = rng.rand(-4f32..4f32)    
+            arr[i] = newRational(rng.rand(-20..20).int64, 20.int64)
-        let mtx = newSMatrixFromArray[5, 5, float32](arr)
+        let mtx = newSMatrixFromArray[5, 5, Rational[int64]](arr)
-        var lu = mtx.clone()        
+        var lu = mtx.clone()
         let pivot = lu.lup()
-        let l = lu.tril(1.0)
+        let l = lu.tril(Rational[int64].one)
         let u = lu.triu()
         let err = pivot * mtx - (l * u)
-        check(err.norm2() < 1e-5)
+        check(err.squared_norm2() == Rational[int64].zero)
 
     test "Linear system solve":
         var rng = initRand(101325)