imported constants in hmatrix and smatrix

- added some meaningful unit tests
- splitted classes into heap-allocated (starting with 'H') and stack-allocated (starting with 'S'), unfortunately there is a lot of duplicated code but I am still unable to find
an elegant solution to use the smae code to deal with both the stack-allocated class and the heap-allocated one

benchmark/benchmark.nim

@@ -0,0 +1,33 @@
+import std/random
+import std/os
+import std/parseutils
+from mmath/rational import Rational, newRational, `$`, zero, abs, `<`, `*`, `-`, `/=`, `+`, `+=`
+from mmath/hmatrix import newHMatrix, clone, lup, lu_solve, `*`, `-`, `$`
+from mmath/hvector import newHvector, `-`, norm
+from mmath/bigint import BigInt, newBigInt, `of`, one, zero, abs, `*`, `mod`, `div`, `-`, `+`, `$`, `==`
+
+let argv = commandLineParams()
+
+let size = block:
+        var n : int = 0
+        if argv.len() > 0: 
+                discard parseInt(argv[0], n) 
+        else: 
+                n = 3
+        n
+
+proc nextInt(r : var Rand, min : int, max : int) : int = min + r.rand(max - min)
+
+var rand = initRand(101325)
+let mtx = block:
+    let valueGenerator = proc(i : int, j : int): Rational[BigInt] =
+        return newRational[BigInt](BigInt.of(rand.nextInt(-1000 * size, 1000 * size).int64), BigInt.of((1000 * size).int64))
+    newHMatrix[Rational[BigInt]](size, size, valueGenerator)
+var lu = mtx.clone()
+let pivot = lu.lup()
+let b = block:
+    let generator = proc(i : int) : Rational[BigInt] = newRational(BigInt.of(rand.nextInt(0, size)), BigInt.of(size))
+    newHVector[Rational[BigInt]](size, generator)
+let x = lu.lu_solve(b, pivot)
+let error = mtx * x - b
+echo $norm(error)

benchmark/config.nims

@@ -0,0 +1 @@
+switch("path", "$projectDir/../src")

matrix.nim

@@ -1,329 +0,0 @@
-from sequtils import newSeqWith, map
-from utils import `...`, `-->`
-from vector import Vector, newVector, createVector
-from options import Option, none
-import future
-from random import randomize, random
-
-type
-  Matrix*[T] = object
-    rows, columns : int 
-    data : seq[T]
-  SMatrix[W, H: static[int], T] =
-    array[1..W, array[1..H, T]]
-  Pivot* = ref object
-    data* : Vector[int]
-    permutations* : int
-  SingularMatrixError* = object of ValueError
-  SizeError* = object of ValueError
-  MatrixRef[T] = ref Matrix[T]
-
-proc `[]`(p : Pivot, index : int) : int = p.data[index]
-proc `[]=`(p : var Pivot, index : int, value : int) = p.data[index] = value
-proc len(p : Pivot) : int = p.data.len
-
-proc newPivot[T](size : int) : Pivot =
-    result = new (Pivot) 
-    result = Pivot(data: newVector[int](size), permutations:0)
-    for i in 0...size:
-        result[i] = i
-
-type AbtractMAtrix = Matrix or SMAtrix
-
-iterator iter_walter[T](m: Matrix[T]) : (int,int, T) {.closure.} =
-  for i in 0...m.rows:
-      for j in 0...m.columns:
-        yield (i, j, m[i,j])
-
-iterator items*[T](m: Matrix[T]): (int,int, T) =
-  for i in 0...m.rows:
-      for j in 0...m.columns:
-        yield (i, j, m[i,j])
-
-proc size*[T](m : Matrix[T]) : (int, int) = (m.rows, m.columns)
-
-proc newMatrix*[T](rows, columns : int, init : T = 0) : Matrix[T] = 
-    result = Matrix[T](rows: rows, columns:columns, data: newSeq[T](rows * columns))
-    for i,j,_ in items(result):
-        result[i,j] = init
-
-proc newMatrix*[T](rows, columns : int, values : openarray[T]) : Matrix[T] = 
-    result = Matrix[T](rows: rows, columns:columns, data: newSeq[T](rows * columns))
-    for i,j,_ in items(result):
-        result[i,j] = values[i * columns + j] 
-
-proc identity*[T](sz : int) : Matrix[T] =
-    result = newMatrix[T](sz,sz,0)
-    for i in 0...sz:
-        result[i,i] = 1
-
-
-proc `[]`*[T](m : Matrix[T], r,c :int) : T =
-    m.data[r*m.columns + c] 
-
-proc `[]`*[T](m : var Matrix[T], r,c :int) : var T =
-    m.data[r*m.columns + c] 
-
-proc `[]=`*[T](m : var Matrix[T], r,c :int, newValue : T) =
-    m.data[r*m.columns + c] = newValue
-
-# proc `[]=`*[T](m : var Matrix[T], r,c :int, newValue : T) =
-#     m.data[r*m.columns + c] = newValue
-
-proc `$`*[T](m : Matrix[T]) : string = 
-    result = "["
-    for i,j,v in items(m):
-        if j == 0:
-            if i > 0: result &= " " 
-            result &= "["
-        result &= $v
-        if j == (m.columns - 1): 
-            result &= "]"
-            if i != (m.rows - 1):
-                result &= ",\n"
-            else:
-                result &= "]\n"
-        else: result &= ", "
-
-proc `*`*[T](m1 : Matrix[T], m2 : Matrix[T]) : Matrix[T] =
-    result = newMatrix[T](m1.rows, m2.columns, 0)
-    for i in 0...result.rows:
-        for j in 0...result.columns:
-            for k in 0...m1.columns:
-                result[i, j] = result[i, j] + m1[i, k] * m2[k, j]
-
-proc `*`*[T](m1 : Matrix[T], v2 : Vector[T]) : Vector[T] =
-    result = newVector[T](m1.rows, 0)
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            result[i] = result[i] + m1[i, j] * v2[j]
-
-proc `+`*[T](m1 : Matrix[T], m2 : Matrix[T]) : Matrix[T] =
-    result = newMatrix[T](m1.rows, m1.columns)
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            result[i,j] = m1[i,j] + m2[i,j]
-
-proc `+`*[T](m1 : Matrix[T], v : T) : Matrix[T] =
-    result = newMatrix[T](m1.rows, m1.columns)
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            result[i,j] = m1[i,j] + v
-
-
-proc `-`*[T](m1 : Matrix[T], m2 : Matrix[T]) : Matrix[T] =
-    result = newMatrix[T](m1.rows, m2.columns)
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            result[i,j] = m1[i,j] - m2[i,j]
-
-proc `-`*[T](m1 : Matrix[T], v : T) : Matrix[T] =
-    result = newMatrix[T](m1.rows, m1.columns)
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            result[i,j] = m1[i,j] - v
-
-proc `+=`*[T](m1 : var Matrix[T], m2 : Matrix[T]) =
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            m1[i,j] += m2[i,j]
-
-proc `+=`*[T](m1 : var Matrix[T], v : T) : Matrix[T] =
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            m1[i,j] += v
-
-proc `-=`*[T](m1 : var Matrix[T], v : T) : Matrix[T] =
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            m1[i,j] -= v
-
-proc `-=`*[T](m1 : var Matrix[T], m2 : Matrix[T]) =
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            m1[i,j] -= m2[i,j]
-
-proc `-`*[T](m : Matrix[T]) : Matrix[T] =
-    result = newMatrix[T](m.rows, m.columns)
-    for i,j,v in m:
-        result[i,j] = -v
-
-proc `==`*[T](m1 : Matrix[T], m2 : Matrix[T]) : bool =
-    if m1.size() != m2.size():
-        return false
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
-            if m1[i,j] != m2[i,j]:
-                return false
-    return true
-
-proc clone*[T](m : Matrix[T]) : Matrix[T] = newMatrix[T](m.rows, m.columns, m.data)
-
-proc transpose*[T](m : Matrix[T]) : Matrix[T] =
-    result = newMatrix[T](m.rows, m.columns)
-    for i, j, v in m:
-        result[j, i] = v
-
-proc det*[T](m : Matrix[T]) : T =
-    var clone = m.clone()
-    clone.gauss_jordan_low()
-    result = 1
-    for i in 0...clone.rows:
-        result *= clone[i, i]
-
-proc swap_rows[T](m : var Matrix[T], id1 : int, id2 : int, pivot : Pivot=nil, other : MatrixRef[T]=nil) =
-    for i in 0...m.columns:
-        let tmp = m[id1, i]
-        m[id1, i] = m[id2, i]
-        m[id2, i] = tmp
-        if other != nil:
-            other[].swap_rows(id1, id2)
-    if pivot != nil:
-        var pv = pivot
-        let tmp = pv[id1]
-        pv[id1] = pv[id2]
-        pv[id2] = tmp
-        pv.permutations =  pv.permutations + 1
-
-proc add_row[T](m : var Matrix[T], sourceIndex : int, destIndex : int, factor : T, other : MatrixRef[T]=nil) =
-    for i in 0...m.columns:
-        m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor
-    if other != nil:
-        other[].add_row(sourceIndex, destIndex, factor)
-
-proc gauss_jordan_low*[T](m : var Matrix[T], other : MatrixRef[T]=nil) =
-    var pivot = newPivot[T](m.rows)
-    for i in 0...m.rows:
-        if m[i, i] == 0:
-            for j in (i + 1)...m.columns:
-                if m[j, i] != 0:
-                    m.swap_rows(i, j, pivot, other)
-                    break
-        for j in (i + 1)...m.rows:
-            if m[i, i] != 0:
-                let factor = -m[j, i] / m[i, i]
-                m.add_row(i, j, factor, other)
-
-proc gauss_jordan_high*[T](m : var Matrix[T], other : MatrixRef[T]=nil) = 
-    var pivot = newPivot[T](m.rows)
-    for i in m.rows-->0:
-        if m[i, i] == 0:
-            for j in i-->0:
-                if m[j, i] != 0:
-                    m.swap_rows(i, j, pivot, other)
-                    break
-        for j in i-->0:
-            if m[i, i] != 0:
-                let factor = -m[j, i] / m[i, i]
-                m.add_row(i, j, factor, other)
-
-proc invert*[T](m : Matrix[T]) : Matrix[T] =
-    var tmp = m.clone()
-    var res = new(MatrixRef[T])
-    res[] = identity[T](tmp.rows)
-    tmp.gauss_jordan_low(res)
-    tmp.gauss_jordan_high(res)
-    for i in 0...res.rows:
-        let f = tmp[i, i]
-        for j in 0...res.columns:
-            res[][i, j] /= f
-    res[]
-
-proc triu*[T](m : Matrix[T], diag_replace=nil) : Matrix[T] =
-    result = Matrix(m.rows, m.columns, 0)
-    for i in range(m.rows):
-        for j in range(i, m.columns):
-            if diag_replace and i == j:
-                result[i, j] = diag_replace
-            else:
-                result[i, j] = m[i, j]
-
-proc tril*[T](m : Matrix[T], diag_replace=nil) : Matrix[T] =
-    result = Matrix(m.rows, m.columns, 0)
-    for i in range(m.rows):
-        for j in range(i + 1):
-            if diag_replace and i == j:
-                result[i, j] = diag_replace
-            else:
-                result[i, j] = m[i, j]
-
-proc lu_row[T](m : var Matrix[T], i : int) =
-    if m[i, i] == 0:
-        raise newException(SingularMatrixError, "Matrix is singular")
-    for j in i...m.columns:
-        for k in 0...i:
-            m[i, j] = m[i, j] - m[i, k] * m[k, j]
-    for j in (i + 1)...m.columns:
-        for k in 0...i:
-            m[j, i] = m[j, i] - m[j, k] * m[k, i]
-        m[j, i] = m[j, i] / m[i, i]
-
-proc lu_pivot[T](m : var Matrix[T], i : int, pivot : Pivot) =
-    var max = abs(m[i, i])
-    var max_index = i
-    for j in (i + 1)...m.rows:
-        if abs(m[j, i]) > max:
-            max = abs(m[i, j])
-            max_index = j
-    if max_index != i:
-        m.swap_rows(i, max_index, pivot)
-
-proc lu*[T](m : var Matrix[T], pivoting=true) : Pivot =
-    var pivot = newPivot[T](m.rows)
-    if pivoting:
-        for i in 0...m.rows:
-            m.lu_pivot(i, pivot)
-            m.lu_row(i)
-    else:
-        for i in 0...m.rows:
-            m.lu_row(i)
-    return pivot
-
-proc lu_solve*[T](m : Matrix[T], b : Vector[T], p : Pivot = nil) : Vector[T] = 
-    var pivot = p
-    if pivot == nil:
-        pivot = new(Pivot)
-        pivot = newPivot[T](m.rows)
-    var x = newVector[T](m.rows)
-    for i in 0...m.rows:
-        x[i] = b[pivot[i]]
-        for k in 0...i:
-            x[i] = x[i] - m[i, k] * x[k]
-    for i in m.rows-->0:
-        for k in (i + 1)...m.rows:
-            x[i] = x[i] - m[i, k] * x[k]
-        x[i] = x[i] / m[i, i]
-    return x
-
-proc lu_invert*[T](m : Matrix[T], pivot : Pivot=nil) : Matrix[T] =
-    if not pivot:
-        pivot = newPivot[T](m.rows)
-    result = newMatrix[T](m.rows, m.columns)
-    for i in 0...m.rows:
-        for j in 0...m.rows:
-            if pivot[j] == i:
-                result[j, i] = 1
-            else:
-                result[j, i] = 0
-            for k in range(j):
-                result[j, i] -= m[j, k] * result[k, i]
-        for j in m.rows-->0:
-            for k in range(j + 1, m.rows):
-                result[j, i] -= m[j, k] * result[k, i]
-            result[j, i] = result[j, i] / m[j, j]
-
-proc lu_det*[T](m : Matrix[T], pivot : Pivot=nil) : T =
-    if not pivot:
-        pivot = newPivot[T](m.rows)
-    result = 1
-    for i in 0...m.rows:
-        result *= m[i, i]
-    if pivot.permutations mod 2 != 0:
-        result *= -1
-
-proc from_pivot[T](pivot : Pivot): Matrix[T] =
-    result = Matrix(len(pivot), len(pivot))
-    for i in 0...pivot.len:
-        result[pivot[i], i] = 1
-        # result[j, i] = 1
-

matrix_test.nim

@@ -1,183 +0,0 @@
-from random import randomize, random
-from utils import `...`
-import matrix
-from matrix import det
-from vector import createVector, newVector, `-`, abs, norm
-import unittest
-
-suite "Nim linear algebra library":
-    echo "suite setup: run once before the tests"
-    randomize()
-    var mtx : Matrix[float64]
-  
-    setup:
-        echo "run before each test"
-        let DIM = 100
-        var numbers = newSeq[float64](DIM * DIM)
-        for i in 0...len(numbers):
-            numbers[i] = float64(random(-DIM..DIM))
-        mtx = newMatrix[float64](DIM,DIM,numbers)
-  
-    teardown:
-        echo "run after each test"
-
-    test "LU decomposition":
-        var lu = mtx.clone()
-        let pivot = lu.lu()
-
-    test "Linear system solve":
-        # let nums = [2.0,1.0,3.0,2.0,6.0,8.0,6.0,8.0,18.0]
-        # let nums = [-3.0, -1.0, 0.0, -2.0, 0.0, 2.0, -3.0, -1.0, 0.0]
-        let DIM = 100
-        var numbers = newSeq[float64](DIM * DIM)
-        for i in 0...len(numbers):
-            numbers[i] = float64(random(-DIM..DIM))
-
-        var mtx = newMatrix[float64](DIM,DIM,numbers)
-        var lu = mtx.clone()
-        let pivot = lu.lu()
-        var b = newVector[float64](DIM)
-        for i in 0...len(b):
-            b[i] = float64(random(DIM))
-
-        let x = lu.lu_solve(b, pivot)
-        let error  = (mtx * x) - b
-
-        check(error.norm() < 1e-5)
-        # var mtx2 = mtx.clone
-        # for i in 0...10000:
-        #     var add  = random(-DIM..DIM).float64()
-        #     discard mtx2 += add
-
-        # echo mtx
-        # mtx[1,1] += 10000.0
-        # echo mtx
-        # echo mtx.det()
-
-        # echo lu
-
-        # give up and stop if this fails
-  
-    echo "suite teardown: run once after the tests"
-
-# let m1 = newMatrix[float](3,3,[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0])
-# var m2 = m1.clone()
-# let m3 = m2.clone()
-# m2[2,1] = -25
-# echo m1 + m2
-# echo m2
-# echo m1.det()
-
-# let nums = [-63, 3, 70, -23, 55,
-#         -100, -37, 81, -98, 84,
-#         -36, -45, -70, 98, -18,
-#         -15, 92, 82, 85, -2,
-#         45, 54, -22, 27, 0
-# ]
-
-# let s = nums.map(proc(n : int) : float = float(n))
-# let s2 = nums.map(n => float(n))
-
-# let m4 = newMatrix[float](5,5,s)
-# echo m4
-# echo m4.det
-
-# var m5 = m4.clone()
-# var pivot = m5.lu()
-# let b = createVector[float](1.0,2.0,3.0,4.0,5.0)
-# let x = m5.lu_solve(b, pivot)
-# echo x
-# proc new(T: typedesc): ref T =
-#   echo "allocating "
-#   new(result)
-
-# var n = new Vector[int]
-
-# type
-#   Index = distinct int
-
-# proc `==` (a, b: Index): bool {.borrow.}
-
-# var af = (0, 0.Index)
-# var b = (0, 0.Index)
-
-# echo af == b # works!
-
-
-# type Person = ref object of RootObj
-#     name : string
-#     age : int
-
-# type Employee = ref object of Person
-#     salary: int
-
-
-# type RecordType = tuple or object
-
-# proc printFields(rec: RecordType) =
-#   for key, value in fieldPairs(rec):
-#     echo key, " = ", value
-
-# proc printFields(rec: ref[RecordType]) =
-#   for key, value in fieldPairs(rec[]):
-#     echo key, " = ", value
-
-
-# let p = Person(name : "Walter", age : 28)
-# let e = Employee(name : "Walter", age : 28, salary: 45000)
-
-# let people : seq[Person] = @[p, e]
-
-# for person in people:
-#     printFields person
-
-# let DIM = 5
-# var numbers = newSeq[float](DIM * DIM)
-# for i in 0...len(numbers):
-#     numbers[i] = float(random(DIM))
-
-# let m = newMatrix[float](DIM,DIM, numbers)
-
-# var m2 = new(Matrix[float])
-# m2[] = m.clone()
-# echo m2[]
-# m2[][0,0] = 300
-# echo m2[]
-
-# let inverse = m2[].invert()
-
-# echo m
-# echo m.det
-# echo m * inverse
-# echo inverse * m
-
-# let nums = [2.0,1.0,3.0,2.0,6.0,8.0,6.0,8.0,18.0]
-# let nums = [-3.0, -1.0, 0.0, -2.0, 0.0, 2.0, -3.0, -1.0, 0.0]
-# let DIM = 600
-# var numbers = newSeq[float64](DIM * DIM)
-# for i in 0...len(numbers):
-#     numbers[i] = float64(random(-DIM..DIM))
-
-# var mtx = newMatrix[float64](DIM,DIM,numbers)
-# var mtx2 = mtx.clone
-# for i in 0...10000:
-#     var add  = random(-DIM..DIM).float64()
-#     discard mtx2 += add
-
-# echo mtx
-# mtx[1,1] += 10000.0
-# echo mtx
-# echo mtx.det()
-
-# var lu = mtx.clone()
-# let pivot = lu.lu()
-# # echo lu
-
-# var b = newVector[float64](DIM)
-# for i in 0...len(b):
-#     b[i] = float64(random(DIM))
-
-# let x = lu.lu_solve(b, pivot)
-# let error  = (mtx * x) - b
-
-# echo abs(error)

mmath.nimble

@@ -0,0 +1,13 @@
+# Package
+
+version       = "0.1.0"
+author        = "Walter Oggioni"
+description   = "Small linear algebra library"
+license       = "MIT"
+srcDir        = "src"
+
+
+# Dependencies
+
+requires "nim >= 0.18"
+requires "nwo >= 0.1"

src/mmath/bigint.nim

@@ -0,0 +1,265 @@
+from error import ParseError
+import nwo/clib
+
+type mpz_t = object
+    mp_alloc : cint
+    mp_size : cint
+    mp_d : pointer 
+
+type mpz_ptr = ptr[mpz_t]
+
+cxface:
+    include "<gmp.h>"
+
+    proc mpz_init(state : mpz_ptr) : void
+    proc mpz_clear(state : mpz_ptr) : void
+    proc mpz_set_si(self : mpz_ptr, value : clong) : void
+    proc mpz_get_str(destination : cstring, base : cint, self : mpz_ptr) : cstring
+    proc mpz_sizeinbase(self : mpz_ptr, base :cint) : csize_t
+    proc mpz_set_str(self : mpz_ptr, std: cstring, base : cint) : int
+    proc mpz_add(rop : mpz_ptr, op1 : mpz_ptr, op2 : mpz_ptr) : void
+    proc mpz_add_ui(rop : mpz_ptr, op1 : mpz_ptr, op2 : culong) : void
+    proc mpz_sub(rop : mpz_ptr, op1 : mpz_ptr, op2 : mpz_ptr) : void
+    proc mpz_sub_ui(rop : mpz_ptr, op : mpz_ptr, op2 : culong) : void
+    proc mpz_ui_sub(rop : mpz_ptr, op2 : culong, op : mpz_ptr) : void
+
+    proc mpz_mul(rop : mpz_ptr, op1 : mpz_ptr, op2 : mpz_ptr) : void
+    proc mpz_mul_si(rop : mpz_ptr, op1 : mpz_ptr, op2 : clong) : void
+    proc mpz_mul_ui(rop : mpz_ptr, op1 : mpz_ptr, op2 : culong) : void
+
+
+    proc mpz_cdiv_q(q : mpz_ptr, n : mpz_ptr, d : mpz_ptr) : void
+    proc mpz_cdiv_r(r : mpz_ptr, n : mpz_ptr, d : mpz_ptr) : void
+    proc mpz_fdiv_q(q : mpz_ptr, n : mpz_ptr, d : mpz_ptr) : void
+    proc mpz_fdiv_r(r : mpz_ptr, n : mpz_ptr, d : mpz_ptr) : void
+    proc mpz_tdiv_q(q : mpz_ptr, n : mpz_ptr, d : mpz_ptr) : void
+    proc mpz_tdiv_r(r : mpz_ptr, n : mpz_ptr, d : mpz_ptr) : void
+
+    proc mpz_tdiv_q_ui(q : mpz_ptr, n : mpz_ptr, d : culong) : void
+    proc mpz_tdiv_r_ui(q : mpz_ptr, n : mpz_ptr, d : culong) : void
+    proc mpz_tdiv_ui(n : mpz_ptr, d : culong) : void
+    proc mpz_neg(rop : mpz_ptr, op : mpz_ptr) : void
+
+    proc mpz_root(rop : mpz_ptr, op1: mpz_ptr, op2: culong) : int
+    proc mpz_sqrt(rop : mpz_ptr, op1: mpz_ptr) : void
+    proc mpz_pow_ui(rop : mpz_ptr, base: mpz_ptr, exp : culong) : void
+    proc mpz_abs(rop : mpz_ptr, op : mpz_ptr) : void
+
+    proc mpz_cmp(op1 : mpz_ptr, op2 : mpz_ptr) : cint
+
+    libs:
+        gmp
+
+proc `=destroy`*(n: mpz_t) =
+    let mpz : mpz_ptr = addr(n)
+    mpz_clear(mpz)
+
+type BigInt* = ref mpz_t
+
+proc newBigInt*(value : int64): BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    mpz_init(mpz)
+    mpz_set_si(mpz, value.clong)
+
+proc toString*(n : BigInt, base : cint) : string =
+    let mpz : mpz_ptr = addr(n[])
+    let size = mpz_sizeinbase(mpz, base.cint)
+    let space : cstring = cast[cstring](alloc(size))
+    discard mpz_get_str(space, base, mpz)
+    result = $space
+    dealloc(space)
+
+proc fromString*(t : type[BigInt], s : string, base : int = 10) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let rc = mpz_set_str(mpz, s, base.cint)
+    if rc != 0:
+        raise newException(ParseError, "Cannot parse integer with base " & $base & " from string " & "'" & s & "'")
+
+proc `$`*(n : BigInt) : string = toString(n, 10)
+
+proc `+`*(n1 : BigInt, n2 : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_init(mpz)
+    mpz_add(mpz, mpz1, mpz2)
+
+proc `+`*(n1 : BigInt, n2 : culong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    mpz_init(mpz)
+    mpz_add_ui(mpz, mpz1, n2.culong)
+
+proc `+`*(n1 : culong, n2 : BigInt) : BigInt = n2 + n1
+
+proc `-`*(n1 : BigInt, n2 : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_init(mpz)
+    mpz_sub(mpz, mpz1, mpz2)
+
+proc `-`*(n1 : BigInt, n2 : culong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    mpz_init(mpz)
+    mpz_sub_ui(mpz, mpz1, n2)
+
+proc `-`*(n1 : culong, n2 : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_init(mpz)
+    mpz_ui_sub(mpz, n1, mpz2)
+
+proc `*`*(n1 : BigInt, n2 : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_init(mpz)
+    mpz_mul(mpz, mpz1, mpz2)
+
+proc `*`*(n1 : BigInt, n2 : clong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    mpz_init(mpz)
+    mpz_mul_si(mpz, mpz1, n2)
+
+proc `*`*(n1 : clong, n2 : BigInt) : BigInt = n2 * n1
+
+proc `*`*(n1 : BigInt, n2 : culong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    mpz_init(mpz)
+    mpz_mul_ui(mpz, mpz1, n2)
+
+proc `*`*(n1 : culong, n2 : BigInt) : BigInt = n2 * n1
+
+proc `div`*(n1 : BigInt, n2 : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_init(mpz)
+    mpz_tdiv_q(mpz, mpz1, mpz2)
+
+proc `div`*(n1 : BigInt, n2 : culong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    mpz_init(mpz)
+    mpz_tdiv_q_ui(mpz, mpz1, n2)
+
+proc `mod`*(n1 : BigInt, n2 : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_init(mpz)
+    mpz_tdiv_r(mpz, mpz1, mpz2)
+
+proc `mod`*(n1 : BigInt, n2 : culong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n1[])
+    mpz_init(mpz)
+    mpz_tdiv_r_ui(mpz, mpz1, n2)
+
+proc `-`*(m : BigInt) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(m[])
+    mpz_neg(mpz, mpz1)
+
+proc `+=`*(n1 : var BigInt, n2 : BigInt) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_add(mpz, mpz, mpz2)
+
+proc `+=`*(n1 : var BigInt, n2 : culong) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    mpz_add_ui(mpz, mpz, n2)
+
+proc `-=`*(n1 : var BigInt, n2 : BigInt) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_sub(mpz, mpz, mpz2)
+
+proc `-=`*(n1 : var BigInt, n2 : culong) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    mpz_sub_ui(mpz, mpz, n2)
+
+proc `*=`*(n1 : var BigInt, n2 : BigInt) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_mul(mpz, mpz, mpz2)
+
+proc `*=`*(n1 : var BigInt, n2 : clong) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    mpz_mul_si(mpz, mpz, n2)
+
+proc `*=`*(n1 : var BigInt, n2 : culong) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    mpz_mul_ui(mpz, mpz, n2)
+
+proc `div=`*(n1 : var BigInt, n2 : BigInt) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_tdiv_q(mpz, mpz, mpz2)
+
+proc `div=`*(n1 : var BigInt, n2 : culong) : void =
+    let mpz : mpz_ptr = addr(n1[])
+    mpz_tdiv_q_ui(mpz, mpz, n2)
+
+proc abs*(n : BigInt) : BigInt = 
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n[])
+    mpz_abs(mpz, mpz1)
+
+proc sqrt*(n : BigInt) : BigInt = 
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n[])
+    mpz_sqrt(mpz, mpz1)
+
+proc pow*(n : BigInt, exp : culong) : BigInt =
+    new(result)
+    let mpz : mpz_ptr = addr(result[])
+    let mpz1 : mpz_ptr = addr(n[])
+    mpz_pow_ui(mpz, mpz1, exp)
+
+proc fact*(n : int) : BigInt =
+    result = newBigInt(1)
+    for i in 1..<n:
+        result *= i
+
+let ZERO = newBigInt(0)
+let ONE = newBigInt(1)
+
+proc one*[BigInt](_: type[BigInt]) : BigInt = ONE
+
+proc zero*[BigInt](_: type[BigInt]) : BigInt = ZERO
+
+proc cmp*(n1 : BigInt, n2 : BigInt) : int = 
+    let mpz1 : mpz_ptr = addr(n1[])
+    let mpz2 : mpz_ptr = addr(n2[])
+    mpz_cmp(mpz1, mpz2).int
+
+proc `<`*(n1 : BigInt, n2 : BigInt) : bool = cmp(n1, n2) < 0
+proc `>`*(n1 : BigInt, n2 : BigInt) : bool = cmp(n1, n2) > 0
+proc `==`*(n1 : BigInt, n2 : BigInt) : bool = cmp(n1, n2) == 0
+proc `<=`*(n1 : BigInt, n2 : BigInt) : bool = cmp(n1, n2) <= 0
+proc `>=`*(n1 : BigInt, n2 : BigInt) : bool = cmp(n1, n2) >= 0
+proc `!=`*(n1 : BigInt, n2 : BigInt) : bool = cmp(n1, n2) != 0
+
+proc `of`*(_ : typedesc[BigInt], n : int64) : BigInt = newBigInt(n)

src/mmath/constant.nim

@@ -0,0 +1,39 @@
+proc one*[T : int64](_: type[T]) : T = 1.T
+
+proc zero*[T : int64](_: type[T]) : T = 0.T
+
+proc one*[T : uint64](_: type[T]) : T = 1.T
+
+proc zero*[T : uint64](_: type[T]) : T = 0.T
+
+proc one*[T : int](_: type[T]) : T = 1.T
+
+proc zero*[T : int](_: type[T]) : T = 0.T
+
+proc one*[T : uint](_: type[T]) : T = 1.T
+
+proc zero*[T : uint](_: type[T]) : T = 0.T
+
+proc one*[T : int16](_: type[T]) : T = 1.T
+
+proc zero*[T : int16](_: type[T]) : T = 0.T
+
+proc one*[T : uint16](_: type[T]) : T = 1.T
+
+proc zero*[T : uint16](_: type[T]) : T = 0.T
+
+proc one*[T : int8](_: type[T]) : T = 1.T
+
+proc zero*[T : int8](_: type[T]) : T = 0.T
+
+proc one*[T : uint8](_: type[T]) : T = 1.T
+
+proc zero*[T : uint8](_: type[T]) : T = 0.T
+
+proc one*[T : float32](_: type[T]) : T = 1.T
+
+proc zero*[T : float32](_: type[T]) : T = 0.T
+
+proc one*[T : float64](_: type[T]) : T = 1.T
+
+proc zero*[T : float64](_: type[T]) : T = 0.T

src/mmath/error.nim

@@ -0,0 +1,4 @@
+type
+    SingularMatrixError* = object of ValueError
+    SizeError* = object of ValueError
+    ParseError* = object of ValueError

src/mmath/hmatrix.nim

@@ -0,0 +1,408 @@
+from nwo/utils import `-->`
+from hvector import HVector, newHVector, buildHVector
+from pivot import HPivot, newHPivot, `[]`, `[]=`, len
+from error import SizeError, SingularMatrixError
+from math import sqrt
+import constant
+
+
+type
+  HMatrix*[T] = object
+    rows, columns : int 
+    data : seq[T]
+
+proc size*[T](m : HMatrix[T]) : (int, int) = (m.rows, m.columns)
+
+proc `[]`*[T](m : HMatrix[T], r, c :int) : T =
+    m.data[r * m.columns + c] 
+
+proc `[]`*[T](m : var HMatrix[T], r, c :int) : var T =
+    m.data[r * m.columns + c] 
+
+proc `[]=`*[T](m : var HMatrix[T], r, c :int, newValue : T) =
+    m.data[r * m.columns + c] = newValue
+
+iterator items*[T](m: HMatrix[T]): (int, int, T) =
+    for i in 0..<m.rows:
+        for j in 0..<m.columns:
+            yield (i, j, m[i, j])
+
+proc rawHMatrix[T](rows, columns : int) : HMatrix[T] = 
+    HMatrix[T](rows: rows, columns:columns, data: newSeq[T](rows * columns))
+
+proc newHMatrix*[T](rows, columns : int, init : proc(i : int, j: int) : T) : HMatrix[T] = 
+    result = rawHMatrix[T](rows, columns)
+    for i in 0..<rows:
+        for j in 0..<columns:
+            result[i,j] = init(i,j)
+
+proc newHMatrix*[T](rows, columns : int, init : T = T.zero) : HMatrix[T] = 
+    result = rawHMatrix[T](rows, columns)
+    for i,j,_ in items(result):
+        result[i,j] = init
+
+proc newHMatrix*[T](rows, columns : int, values : openarray[T]) : HMatrix[T] = 
+    result = rawHMatrix[T](rows, columns)
+    for i, j, _ in items(result):
+        result[i,j] = values[i * columns + j] 
+
+proc identity*[T](sz : int) : HMatrix[T] =
+    let init = proc(i : int, j: int) : T =
+        if i == j:
+            result = T.one
+        else:
+            result = T.zero
+    result = newHMatrix[T](sz, sz, init)
+
+proc `$`*[T](m : HMatrix[T]) : string = 
+    result = "["
+    for i,j,v in items(m):
+        if j == 0:
+            if i > 0: result &= " " 
+            result &= "["
+        result &= $v
+        if j == (m.columns - 1): 
+            result &= "]"
+            if i != (m.rows - 1):
+                result &= ",\n"
+            else:
+                result &= "]\n"
+        else: result &= ", "
+
+proc `*`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] =
+    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:
+                result[i, j] = result[i, j] + m1[i, k] * m2[k, j]
+
+proc `*`*[T](m1 : HMatrix[T], v2 : HVector[T]) : HVector[T] =
+    result = newHVector[T](m1.rows)
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i] = result[i] + m1[i, j] * v2[j]
+
+proc `+`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] =
+    result = newHMatrix[T](m1.rows, m1.columns)
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] + m2[i,j]
+
+proc `+`*[T](m1 : HMatrix[T], v : T) : HMatrix[T] =
+    result = newHMatrix[T](m1.rows, m1.columns)
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] + v
+
+
+proc `-`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] =
+    result = newHMatrix[T](m1.rows, m2.columns)
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] - m2[i,j]
+
+proc `-`*[T](m1 : HMatrix[T], v : T) : HMatrix[T] =
+    result = newHMatrix[T](m1.rows, m1.columns)
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] - v
+
+proc `+=`*[T](m1 : var HMatrix[T], m2 : HMatrix[T]) =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] += m2[i,j]
+
+proc `+=`*[T](m1 : var HMatrix[T], v : T) : HMatrix[T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] += v
+
+proc `-=`*[T](m1 : var HMatrix[T], v : T) : HMatrix[T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] -= v
+
+proc `-=`*[T](m1 : var HMatrix[T], m2 : HMatrix[T]) =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] -= m2[i,j]
+
+proc `-`*[T](m : HMatrix[T]) : HMatrix[T] =
+    result = newHMatrix[T](m.rows, m.columns)
+    for i,j,v in m:
+        result[i,j] = -v
+
+proc `==`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : bool =
+    if m1.size() != m2.size():
+        return false
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            if m1[i,j] != m2[i,j]:
+                return false
+    return true
+
+proc clone*[T](m : HMatrix[T]) : HMatrix[T] = newHMatrix[T](m.rows, m.columns, m.data)
+
+proc transpose*[T](m : HMatrix[T]) : HMatrix[T] =
+    result = newHMatrix[T](m.columns, m.rows)
+    for i, j, v in items(m):
+        result[j, i] = v
+
+proc swap_rows[T](m : var HMatrix[T], id1 : int, id2 : int) =
+    for i in 0..<m.columns:
+        let tmp = m[id1, i]
+        m[id1, i] = m[id2, i]
+        m[id2, i] = tmp
+
+proc swap_rows[T](m : var HMatrix[T], id1 : int, id2 : int, pivot : var HPivot[T]) =
+    m.swap_rows(id1, id2)
+    let tmp = pivot[id1]
+    pivot[id1] = pivot[id2]
+    pivot[id2] = tmp
+    pivot.permutations += 1
+        
+proc swap_rows[T](
+        m : var HMatrix[T], id1 : int, 
+        id2 : int, 
+        pivot : var HPivot[T], 
+        other : var HMatrix[T]) =
+    m.swap_rows(id1, id2, pivot)
+    other.swap_rows(id1, id2)
+
+proc add_row[T](
+        m : var HMatrix[T], 
+        sourceIndex : int, 
+        destIndex : int, 
+        factor : T) =
+    for i in 0..<m.columns:
+        m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor
+
+proc add_row[T](
+        m : var HMatrix[T], 
+        sourceIndex : int, 
+        destIndex : int, 
+        factor : T, 
+        other : var HMatrix[T]) =
+    add_row(m, source_index, dest_index, factor)
+    other.add_row(sourceIndex, destIndex, factor)
+
+proc gauss_jordan_low*[T](
+        m : var HMatrix[T],
+        other : var HMatrix[T]) =
+    var pivot = newHPivot[T](m.rows)
+    for i in 0..<m.rows:
+        if m[i, i] == 0:
+            for j in (i + 1)..<m.columns:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot, other)
+                    break
+        for j in (i + 1)..<m.rows:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor, other)
+
+proc gauss_jordan_low*[T](
+        m : var HMatrix[T]) =
+    var pivot = newHPivot[T](m.rows)
+    for i in 0..<m.rows:
+        if m[i, i] == 0:
+            for j in (i + 1)..<m.columns:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot)
+                    break
+        for j in (i + 1)..<m.rows:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor)
+
+proc gauss_jordan_high*[T](
+        m : var HMatrix[T]) =
+    var pivot = newHPivot[T]()
+    for i in m.rows-->0:
+        if m[i, i] == 0:
+            for j in i-->0:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot)
+                    break
+        for j in i-->0:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor)
+
+proc gauss_jordan_high*[T](
+        m : var HMatrix[T],
+        other : var HMatrix[T]) =
+    var pivot = newHPivot[T](m.rows)
+    for i in m.rows-->0:
+        if m[i, i] == 0:
+            for j in i-->0:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot, other)
+                    break
+        for j in i-->0:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor, other)
+
+proc det*[T](m : HMatrix[T]) : T =
+    if m.rows != m.columns:
+        raise newException(SizeError, "Matrix must be square in order to compute the determinant")
+    var clone = m.clone()
+    clone.gauss_jordan_low()
+    result = T.one
+    for i in 0..<clone.rows:
+        result *= clone[i, i]
+
+proc invert*[T](m : HMatrix[T]) : HMatrix[T] =
+    if m.rows != m.columns:
+        raise newException(SizeError, "Matrix must be square in order to compute the determinant")
+    var tmp = m.clone()
+    result = identity[T](m.rows)
+    tmp.gauss_jordan_low(result)
+    tmp.gauss_jordan_high(result)
+    for i in 0..<result.rows:
+        let f = tmp[i, i]
+        for j in 0..<result.columns:
+            result[i, j] /= f
+
+proc triu*[T](m : HMatrix[T], diag_replace: T) : HMatrix[T] =
+    result = newHMatrix[T](m.rows, m.columns, T.zero)
+    for i in 0..<m.rows:
+        for j in i..<m.columns:
+            if i == j:
+                result[i, j] = diag_replace
+            else:
+                result[i, j] = m[i, j]
+
+proc triu*[T](m : HMatrix[T]) : HMatrix[T] =
+    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, T.zero)
+    for i in 0..<m.rows:
+        for j in 0..<(i + 1):
+            if i == j:
+                result[i, j] = diag_replacement
+            else:
+                result[i, j] = m[i, j]
+
+proc tril*[T](m : HMatrix[T]) : HMatrix[T] =
+    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] == T.zero:
+        raise newException(SingularMatrixError, "Matrix is singular")
+    for j in i..<m.columns:
+        for k in 0..<i:
+            m[i, j] = m[i, j] - m[i, k] * m[k, j]
+    for j in (i + 1)..<m.columns:
+        for k in 0..<i:
+            m[j, i] = m[j, i] - m[j, k] * m[k, i]
+        m[j, i] /= m[i, i]
+
+proc lu_pivot[T](m : var HMatrix[T], i : int, pivot : var HPivot[T]) =
+    var max = abs(m[i, i])
+    var max_index = i
+    for j in (i + 1)..<m.rows:
+        if abs(m[j, i]) > max:
+            max = abs(m[i, j])
+            max_index = j
+    if max_index != i:
+        m.swap_rows(i, max_index, pivot)
+
+proc lup*[T](m : var HMatrix[T]) : HPivot[T] =
+    result = newHPivot[T](m.rows)
+    for i in 0..<m.rows:
+        m.lu_pivot(i, result)
+        m.lu_row(i)
+
+proc lu*[T](m : var HMatrix[T]) =
+    for i in 0..<m.rows:
+        m.lu_row(i)
+
+proc lu_solve*[T](m : HMatrix[T], b : HVector[T], pivot : HPivot[T]) : HVector[T] = 
+    var x = newHVector[T](m.rows)
+    for i in 0..<m.rows:
+        x[i] = b[pivot[i]]
+        for k in 0..<i:
+            x[i] = x[i] - m[i, k] * x[k]
+    for i in m.rows --> 0:
+        for k in (i + 1)..<m.rows:
+            x[i] = x[i] - m[i, k] * x[k]
+        if m[i,i] != T.zero:
+            x[i] /= m[i, i]
+        else:
+            raise newException(SingularMatrixError, "Matrix is singular")
+    return x
+
+proc lu_solve*[T](
+        m : HMatrix[T], 
+        b : HVector[T]) : HVector[T] = 
+    var pivot = newHPivot[T]()
+    lu_solve(m, b, pivot)
+
+proc lu_invert*[T](m : HMatrix[T], pivot : HPivot[T]) : HMatrix[T] =
+    if m.rows != m.columns:
+        raise newException(SizeError, "Matrix must be square in order to compute the inverse")
+    result = newHMatrix[T](m.rows, m.columns)
+    for i in 0..<m.rows:
+        for j in 0..<m.rows:
+            if pivot[j] == i:
+                result[j, i] = T.one
+            else:
+                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:
+            for k in range(j + 1, m.rows):
+                result[j, i] -= m[j, k] * result[k, i]
+            result[j, i] = result[j, i] / m[j, j]
+
+proc lu_invert*[T](m : HMatrix[T]) : HMatrix[T] = lu_invert(m, newHPivot[T]())
+
+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 = T.one
+    for i in 0..<m.rows:
+        result *= m[i, i]
+    if pivot.permutations mod 2 != 0:
+        result *= -T.one
+
+proc lu_det*[T](m : HMatrix[T]) : T = 
+    if m.rows != m.columns:
+        raise newException(SizeError, "Matrix must be square in order to compute the determinant")
+    var clone = m.clone()
+    lu_det(clone)
+
+proc squared_norm2*[T](m : HMatrix[T]): T =
+    result = T.zero
+    for i, j, v in items(m):
+        result += v * v
+
+proc norm2*[T](m : HMatrix[T]): T =
+    sqrt(m.squared_norm2())
+
+proc `*`*[T](pivot : HPivot[T], m : HMatrix[T]) : HMatrix[T] =
+    result = m.clone()
+    var pclone = pivot
+    for i in 0..<pclone.len():
+        while i != pclone[i]:
+            result.swap_rows(i, pclone[i])
+            let tmp = pclone[i]
+            pclone[i] = pclone[tmp]
+            pclone[tmp] = tmp
+
+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] = T.one
+

src/mmath/hvector.nim

@@ -0,0 +1,50 @@
+type HVector*[T] = seq[T]
+
+proc newHVector*[T](size : int, init: T=T.zero) : HVector[T] = 
+    result = newSeq[T](size)
+    for i in 0..<len(result):
+        result[i] = init
+
+proc newHVector*[T](size : int, init: proc(index : int) : T) : HVector[T] = 
+    result = newSeq[T](size)
+    for i in 0..<len(result):
+        result[i] = init(i)
+
+proc `+`*[T](v1 : HVector[T], v2:HVector[T]) : HVector[T] =
+    result = newHVector[T](len(v1))
+    for i in 0..<len(v1):
+        result[i] = v1[i] + v2[i]
+
+proc `-`*[T](v1 : HVector[T], v2:HVector[T]) : HVector[T] =
+    result = newHVector[T](len(v1))
+    for i in 0..<len(v1):
+        result[i] = v1[i] - v2[i]
+
+proc `*`*[T](v1 : HVector[T], v2:HVector[T]) : T =
+    result = T.zero
+    for i in 0..<len(v1):
+        result += v1[i] * v2[i]
+
+proc `+=`*[T](v1 : var HVector[T], v2 : HVector[T]) = 
+    for i in 0..<len(v1):
+        v1[i] += v2[i]
+
+proc `-=`*[T](v1 : var HVector[T], v2 : HVector[T]) = 
+    for i in 0..<len(v1):
+        v1[i] -= v2[i]
+
+proc `+=`*[T](v : var HVector[T], value : T) = v.add(value)
+
+proc buildHVector*[T](elems : varargs[T]) : HVector[T] = 
+    result = newSeq[T]()
+    for elem in items(elems):
+        result += elem
+
+proc norm*[T](v : HVector[T]) : T =
+    result = T.zero
+    for value in v:
+        result += v * v
+
+proc abs*[T](v : HVector[T]) : T =
+    return v.norm().sqrt()
+

src/mmath/pivot.nim

@@ -0,0 +1,27 @@
+from hvector import HVector, newHVector
+
+type
+    HPivot*[T] = object
+      data : HVector[int]
+      permutations* : int
+
+    SPivot*[SIZE : static[int], T] = object
+      data : array[SIZE, int]
+      permutations* : int
+  
+proc `[]`*[T](p : HPivot[T], index : int) : int = p.data[index]
+proc `[]=`*[T](p : var HPivot[T], index : int, value : int) = p.data[index] = value
+proc len*[T](p : HPivot[T]) : int = p.data.len
+proc `$`*[T](pivot : HPivot[T]) : string = $pivot.data
+proc newHPivot*[T](size : int) : HPivot[T] =
+    result = HPivot[T](data: newHVector[int](size), permutations:0)
+    for i in 0..<size:
+        result[i] = i
+
+proc `[]`*[SIZE, T](p : SPivot[SIZE, T], index : int) : int = p.data[index]
+proc `[]=`*[SIZE, T](p : var SPivot[SIZE, T], index : int, value : int) = p.data[index] = value
+proc len*[SIZE, T](p : SPivot[SIZE, T]) : int = SIZE
+proc `$`*[SIZE, T](pivot : SPivot[SIZE, T]) : string = $pivot.data
+proc newSPivot*[SIZE, T]() : SPivot[SIZE, T] =
+    for i in 0..<SIZE:
+        result[i] = i

src/mmath/rational.nim

@@ -0,0 +1,118 @@
+import constant
+
+export constant
+
+proc gcd*[T](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](n1 : T, n2 : T) : T = n1 * n2 div gcd(n1, n2)
+
+proc pow*[T](n : T, e : T) : T =
+    result = T.one
+    for i in 0..<e:
+        result *= n    
+
+type Rational*[T] = 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.one)
+proc newRational*[T](num : T, den: T) : Rational[T] = Rational[T](num: num, den: den)
+
+proc zero*[T](_: type[Rational[T]]): Rational[T] = newRational(T.zero, T.one)
+
+proc one*[T](_: type[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

@@ -0,0 +1,391 @@
+from nwo/utils import `-->`, box
+from svector import SVector
+from pivot import SPivot, newSPivot, `[]`, `[]=`, len
+from error import SizeError, SingularMatrixError
+from math import sqrt
+import constant
+
+type
+  SMatrix*[ROWS, COLUMNS: static[int], T] = object
+    data : array[0..(ROWS*COLUMNS - 1), T]
+  SquareSMatrix*[SIZE: static[int], T] = SMatrix[SIZE, SIZE, T]
+
+proc size*[ROWS, COLUMNS : static[int], T](m : SMatrix) : (int, int) = (m.rows, m.columns)
+
+proc rows*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : int = ROWS
+proc columns*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : int = COLUMNS
+
+proc `[]`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS,COLUMNS,T], r, c :int) : T = 
+    m.data[r * COLUMNS + c]
+
+proc `[]`*[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS,COLUMNS,T], r, c :int) : var T = 
+    m.data[r * COLUMNS + c]
+
+proc `[]=`*[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS,COLUMNS,T], r, c :int, newValue : T) =
+    m.data[r * COLUMNS + c] = newValue
+
+iterator items*[ROWS, COLUMNS : static[int], T](m: SMatrix[ROWS,COLUMNS,T]): (int, int, T) =
+    for i in 0..<m.rows:
+        for j in 0..<m.columns:
+            yield (i, j, m[i, j])
+
+proc newSMatrix*[ROWS, COLUMNS : static[int], T](init : T) : SMatrix[ROWS, COLUMNS, T] =
+    for i,j,_ in items(result):
+        result[i,j] = init
+
+proc newSMatrix*[ROWS, COLUMNS : static[int], T](init : proc(i : int, j : int)) : SMatrix[ROWS, COLUMNS, T] =
+    for i in 0..<ROWS:
+        for j in 0..<COLUMNS:
+            result[i,j] = init(i, j)
+
+proc newSMatrixFromArray*[ROWS, COLUMNS : static[int], T](values : array[0..(ROWS * COLUMNS - 1), T]) : auto = 
+    SMatrix[ROWS, COLUMNS, T](data:values)
+
+proc identity*[SIZE: static[int], T]() : SquareSMatrix[SIZE, T] =
+    for i in 0..<SIZE:
+        result[i,i] = T.one
+
+proc `$`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : string = 
+    result = "["
+    for i,j,v in items(m):
+        if j == 0:
+            if i > 0: result &= " " 
+            result &= "["
+        result &= $v
+        if j == (m.columns - 1): 
+            result &= "]"
+            if i != (m.rows - 1):
+                result &= ",\n"
+            else:
+                result &= "]\n"
+        else: result &= ", "
+
+proc `*`*[ROWS1, COLUMNS2, COMMON : static[int], T](
+        m1 : SMatrix[ROWS1, COMMON, 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] += 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:
+        for j in 0..<m.columns:
+            result[j] += m[i, j] * v[i]
+
+proc `*`*[SIZE : static[int], T](m1 : SquareSMatrix[SIZE, T], v2 : SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i] += m1[i, j] * v2[j]
+
+proc `+`*[ROWS, COLUMNS : static[int], T](addend1 : SMatrix[ROWS, COLUMNS, T], addend2 : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] =
+    for i in 0..<addend1.rows:
+        for j in 0..<addend1.columns:
+            result[i,j] = addend1[i,j] + addend2[i,j]
+
+proc `+`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] + v
+
+proc `-`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS,T], m2 : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS,T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] - m2[i,j]
+
+proc `-`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            result[i,j] = m1[i,j] - v
+
+proc `+=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], m2 : SMatrix[ROWS, COLUMNS, T]) =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] += m2[i,j]
+
+proc `+=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] += v
+
+proc `-=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] -= v
+
+proc `-=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], m2 : SMatrix[ROWS, COLUMNS, T]) =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            m1[i,j] -= m2[i,j]
+
+proc `-`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] =
+    for i,j,v in m:
+        result[i,j] = -v
+
+proc `==`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], m2 : SMatrix[ROWS, COLUMNS, T]) : bool =
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
+            if m1[i,j] != m2[i,j]:
+                return false
+    return true
+
+proc squared_norm2*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]): T =
+    result = T.zero
+    for i, j, v in items(m):
+        result += v * v
+
+proc norm2*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]): T =
+    sqrt(m.squared_norm2())
+
+proc clone*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] = 
+    newSMatrixFromArray[ROWS, COLUMNS, T](m.data)
+
+proc transpose*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[COLUMNS, ROWS, T] =
+    for i, j, v in items(m):
+        result[j, i] = v
+
+proc swap_rows[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS, COLUMNS, T], id1 : int, id2 : int) =
+    for i in 0..<m.columns:
+        let tmp = m[id1, i]
+        m[id1, i] = m[id2, i]
+        m[id2, i] = tmp
+
+proc swap_rows[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS, COLUMNS, T], id1 : int, id2 : int, pivot : var SPivot[ROWS, T]) =
+    m.swap_rows(id1, id2)
+    let tmp = pivot[id1]
+    pivot[id1] = pivot[id2]
+    pivot[id2] = tmp
+    pivot.permutations += 1
+        
+proc swap_rows[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T], id1 : int, 
+        id2 : int, 
+        pivot : var SPivot[ROWS, T], 
+        other : var SMatrix[ROWS, COLUMNS,T]) =
+    m.swap_rows(id1, id2, pivot)
+    other.swap_rows(id1, id2)
+
+proc add_row[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T], 
+        sourceIndex : int, 
+        destIndex : int, 
+        factor : T) =
+    for i in 0..<m.columns:
+        m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor
+
+proc add_row[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T], 
+        sourceIndex : int, 
+        destIndex : int, 
+        factor : T, 
+        other : var SMatrix[ROWS, COLUMNS, T]) =
+    add_row(m, source_index, dest_index, factor)
+    other.add_row(sourceIndex, destIndex, factor)
+
+proc gauss_jordan_low*[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T],
+        other : var SMatrix[ROWS, COLUMNS, T]) =
+    var pivot = newSPivot[ROWS, T]()
+    for i in 0..<m.rows:
+        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)
+                    break
+        for j in (i + 1)..<m.rows:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor, other)
+
+proc gauss_jordan_low*[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T]) =
+    var pivot = newSPivot[ROWS, T]()
+    for i in 0..<m.rows:
+        if m[i, i] == 0:
+            for j in (i + 1)..<m.columns:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot)
+                    break
+        for j in (i + 1)..<m.rows:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor)
+
+proc gauss_jordan_high*[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T]) =
+    var pivot = newSPivot[ROWS, T]()
+    for i in m.rows-->0:
+        if m[i, i] == 0:
+            for j in i-->0:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot)
+                    break
+        for j in i-->0:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor)
+
+proc gauss_jordan_high*[ROWS, COLUMNS : static[int], T](
+        m : var SMatrix[ROWS, COLUMNS, T],
+        other : var SMatrix[ROWS, COLUMNS, T]) =
+    var pivot = newSPivot[ROWS, T]()
+    for i in m.rows-->0:
+        if m[i, i] == 0:
+            for j in i-->0:
+                if m[j, i] != 0:
+                    m.swap_rows(i, j, pivot, other)
+                    break
+        for j in i-->0:
+            if m[i, i] != 0:
+                let factor = -m[j, i] / m[i, i]
+                m.add_row(i, j, factor, other)
+
+proc det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T =
+    var clone = m.clone()
+    clone.gauss_jordan_low()
+    result = T.one
+    for i in 0..<clone.rows:
+        result *= clone[i, i]
+
+proc invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] =
+    var tmp = m.clone()
+    result = identity[SIZE, T]()
+    tmp.gauss_jordan_low(result)
+    tmp.gauss_jordan_high(result)
+    for i in 0..<result.rows:
+        let f = tmp[i, i]
+        for j in 0..<result.columns:
+            result[i, j] /= f
+
+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 0..<m.columns:
+            if i < j:
+                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 0..<m.columns:
+            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..<m.columns:
+            if i > j:
+                result[i, j] = m[i, j]
+            elif i == j:
+                result[i, j] = diag_replacement
+            else:
+                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..<m.columns:
+            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] == T.zero:
+        raise newException(SingularMatrixError, "Matrix is singular")
+    for j in i..<m.columns:
+        for k in 0..<i:
+            m[i, j] = m[i, j] - m[i, k] * m[k, j]
+    for j in (i + 1)..<m.columns:
+        for k in 0..<i:
+            m[j, i] = m[j, i] - m[j, k] * m[k, i]
+        m[j, i] /= m[i, i]
+
+proc lu_pivot[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T], i : int, pivot : var SPivot[SIZE, T]) =
+    var max = abs(m[i, i])
+    var max_index = i
+    for j in (i + 1)..<m.rows:
+        if abs(m[j, i]) > max:
+            max = abs(m[i, j])
+            max_index = j
+    if max_index != i:
+        m.swap_rows(i, max_index, pivot)
+
+proc lu*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) =
+    for i in 0..<m.rows:
+        m.lu_row(i)
+
+proc lup*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) : SPivot[SIZE, T] =
+    result = newSPivot[SIZE,T]()
+    for i in 0..<m.rows:
+        m.lu_pivot(i, result)
+        m.lu_row(i)            
+
+proc lu_solve*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], b : SVector[SIZE, T], pivot : SPivot[SIZE, T]) : SVector[SIZE, T] = 
+    var x : SVector[SIZE, T]
+    for i in 0..<m.rows:
+        x[i] = b[pivot[i]]
+        for k in 0..<i:
+            x[i] = x[i] - m[i, k] * x[k]
+    for i in m.rows --> 0:
+        for k in (i + 1)..<m.rows:
+            x[i] = x[i] - m[i, k] * x[k]
+        if m[i,i] != 0:
+            x[i] /= m[i, i]
+        else:
+            raise newException(SingularMatrixError, "Matrix is singular")
+    return x
+
+proc lu_solve*[SIZE : static[int], T](
+        m : SquareSMatrix[SIZE, T], 
+        b : SVector[SIZE, T]) : SVector[SIZE, T] = 
+    var pivot = newSPivot[SIZE, T]()
+    lu_solve(m, b, pivot)
+
+proc lu_invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], pivot : SPivot[SIZE, T]) : SquareSMatrix[SIZE, T] =
+    for i in 0..<m.rows:
+        for j in 0..<m.rows:
+            if pivot[j] == i:
+                result[j, i] = T.one
+            else:
+                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:
+            for k in range(j + 1, m.rows):
+                result[j, i] -= m[j, k] * result[k, i]
+            result[j, i] = result[j, i] / m[j, j]
+
+proc lu_invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] = lu_invert(m, newSPivot[SIZE, T]())
+
+proc lu_det*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) : T =
+    let pivot = m.lup()
+    result = T.one
+    for i in 0..<m.rows:
+        result *= m[i, i]
+    if pivot.permutations mod 2 != 0:
+        result *= -T.one
+
+proc lu_det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T = 
+    var clone = m.clone()
+    lu_det(clone)
+
+proc `*`*[ROWS, COLUMNS : static[int], T](pivot : SPivot[ROWS, T], m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] =
+    result = m.clone()
+    var pclone = pivot
+    for i in 0..<pclone.len():
+        while i != pclone[i]:
+            result.swap_rows(i, pclone[i])
+            let tmp = pclone[i]
+            pclone[i] = pclone[tmp]
+            pclone[tmp] = tmp
+
+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] = T.one
+

src/mmath/svector.nim

@@ -0,0 +1,79 @@
+from math import sqrt
+
+type SVector*[S : static[int], T] = array[0..(S-1), T]
+
+proc newSVector*[SIZE, T](init: T=T.zero) : SVector[SIZE, T] = 
+    for i in 0..<len(result):
+        result[i] = init
+
+proc newSVector*[SIZE, T](init: proc(index : int) : T) : SVector[SIZE, T] = 
+    for i in 0..<len(result):
+        result[i] = init(i)
+
+proc buildSVector*[SIZE, T](elems : varargs[T]) : SVector[SIZE, T] = 
+    for i in 0..<elems.len:
+        result[i] = elems[i]
+
+proc `+`*[SIZE, T](v1 : SVector[SIZE, T], scalar: T) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = v1[i] + scalar
+
+proc `+`*[SIZE, T](scalar: T, v1 : SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = scalar + v1[i]        
+        
+proc `-`*[SIZE, T](v1 : SVector[SIZE, T], scalar: T) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = v1[i] - scalar
+
+proc `-`*[SIZE, T](scalar: T, v1 : SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = scalar - v1[i]        
+
+proc `*`*[SIZE, T](v1 : SVector[SIZE, T], scalar: T) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = v1[i] * scalar
+
+proc `*`*[SIZE, T](scalar: T, v1 : SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = scalar * v1[i]
+        
+proc `/`*[SIZE, T](v1 : SVector[SIZE, T], scalar: T) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = v1[i] / scalar
+
+proc `/`*[SIZE, T](scalar: T, v1 : SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = scalar / v1[i]        
+
+proc `+`*[SIZE, T](v1 : SVector[SIZE, T], v2: SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = v1[i] + v2[i]
+
+proc `-`*[SIZE, T](v1 : SVector[SIZE, T], v2: SVector[SIZE, T]) : SVector[SIZE, T] =
+    for i in 0..<len(v1):
+        result[i] = v1[i] - v2[i]
+
+proc `*`*[SIZE, T](v1 : SVector[SIZE, T], v2: SVector[SIZE, T]) : T =
+    result = T.zero
+    for i in 0..<len(v1):
+        result += v1[i] * v2[i]
+
+proc `+=`*[SIZE, T](v1 : var SVector[SIZE, T], v2: SVector[SIZE, T]) =
+    for i in 0..<len(v1):
+        v1[i] += v2[i]
+
+proc `-=`*[SIZE, T](v1 : var SVector[SIZE, T], v2: SVector[SIZE, T]) =
+    for i in 0..<len(v1):
+        v1[i] -= v2[i]
+
+proc `+=`*[SIZE, T](v : var SVector[SIZE, T], value : T) = v.add(value)
+
+proc norm*[SIZE, T](v : SVector[SIZE, T]) : T =
+    result = T.zero
+    for value in v:
+        result += v * v
+
+proc abs*[SIZE, T](v : SVector[SIZE, T]) : T =
+    return math.sqrt(v.norm)
+

tests/config.nims

@@ -0,0 +1 @@
+switch("path", "$projectDir/../src")

tests/test_bigint.nim

@@ -0,0 +1,63 @@
+import mmath/bigint
+import unittest
+import random
+from strutils import parseBiggestUInt
+
+suite "Nim arbitrary precision integers (powered by GMP)":
+
+    let str = "2347822319"
+    test "parse " & str:
+        let str = "2347822319"
+        let bi = BigInt.fromString(str, 10)
+        let ul : BiggestUInt = parseBiggestUInt(str)
+        check newBigInt(ul.int64) == bi
+
+    var rng = initRand(101325)
+
+    for i in 0..5:
+        let n1 = rng.rand(-1000..1000).int64
+        let n2 = rng.rand(-1000..1000).int64
+        let bi1 = BigInt.of(n1)
+        let bi2 = BigInt.of(n2)
+
+        test $n1 & " + " & $n2:
+            check BigInt.of(n1 + n2) == bi1 + bi2
+
+        test $n1 & " - " & $n2:
+            check BigInt.of(n1 - n2) == bi1 - bi2
+
+        test $n1 & " * " & $n2:
+            check BigInt.of(n1 * n2) == bi1 * bi2
+
+        test $n1 & " div " & $n2:
+            check BigInt.of(n1 div n2) == bi1 div bi2
+
+        test $n1 & " mod " & $n2:
+            check BigInt.of(n1 mod n2) == bi1 mod bi2
+        
+        test $n1 & " pow " & $3:
+            check BigInt.of(n1 * n1 * n1) == pow(bi1, 3)
+
+        test "abs(" & $n1 & ")":
+            check BigInt.of(n1.abs()) == bi1.abs() 
+
+        test "abs(" & $n2 & ")":
+            check BigInt.of(n2.abs()) == bi2.abs() 
+
+        test $n1 & " cmp " & $n2:
+            check cmp(n1, n2) == cmp(bi1, bi2)
+
+        test $n1 & " < " & $n2:
+            check (n1 < n2) == (bi1 < bi2)
+
+        test $n1 & " > " & $n2:
+            check (n1 > n2) == (bi1 > bi2)
+
+        test $n1 & " == " & $n1:
+            check bi1 == bi1
+
+        test $n2 & " == " & $n2:
+            check bi2 == bi2
+
+        test $n1 & " != " & $n2:
+            check bi1 != bi2

tests/test_hmatrix.nim

@@ -0,0 +1,115 @@
+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, squared_norm2
+from mmath/hvector import buildHVector, newHVector, `-`, abs, norm
+import unittest
+
+suite "Nim linear algebra library":
+    let test_matrix = newHMatrix[float32](3, 3, [1.0f32, 2f32, 3f32, 4f32, 5f32, 6f32, 8f32, 7f32, 9f32])
+
+    test "+":
+        let mtx1 = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
+        let mtx2 = newHMatrix[int8](3, 3, [-1i8,-2,-3,-4,-5,-6,-7,-8,-9])
+        check(mtx1 + mtx2 == newHMatrix[int8](3,3))
+
+    test "-":
+        let mtx  = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
+        check(mtx - mtx == newHMatrix[int8](3, 3))
+
+    test "+=":
+        var mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
+        mtx += mtx
+        check(mtx == newHMatrix[int8](3, 3, [2i8,4,6,8,10,12,14,16,18]))
+
+    test "-=":
+        var mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
+        mtx -= mtx
+        check(mtx == newHMatrix[int8](3, 3))
+
+    test "*":
+        block:
+            let mtx  = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,8,7,9])
+            let res = mtx * buildHVector[int8](1i8, 2, 3)
+            check(res == buildHVector[int8](14i8, 32, 49))
+        block:
+            let mtx = newHMatrix[int](2, 5, [1,4,8,2,5,7,3,6,9,0])
+            let res = mtx * mtx.transpose()
+            check(res == newHMatrix[int](2, 2, [110,85,85,175]))
+        block:
+            let mtx = newHMatrix[int](2, 5, [1,4,8,2,5,7,3,6,9,0])
+            let res = mtx.transpose() * mtx
+            check(res == newHMatrix[int](5, 5, [
+                50, 25, 50, 65, 5,
+                25, 25, 50, 35, 20,
+                50, 50, 100, 70, 40,
+                65, 35, 70, 85, 10,
+                5, 20, 40, 10, 25]))
+        
+    test "Determinant":
+        check(test_matrix.det() == -9f32)
+
+    test "LU Determinant":
+        check(test_matrix.lu_det() == -9f32)
+
+    test "Inverse":
+        let err = test_matrix * test_matrix.invert() - identity[float32](3)
+        check(err.norm2() < 1e-5)
+
+    test "LU decomposition":
+        var rng = initRand(101325)
+        var arr : array[0..(25 - 1), Rational[int64]]
+        for i in 0..<arr.len:
+            arr[i] = newRational(rng.rand(-20..20).int64, 20.int64)
+        let mtx = newHMatrix[Rational[int64]](5, 5, arr)
+        var lu = mtx.clone()        
+        let pivot = lu.lup()
+        let l = lu.tril(Rational[int64].one)
+        let u = lu.triu()
+        let err = pivot * mtx - (l * u)
+        check(err.squared_norm2() == Rational[int64].zero)
+
+    test "Linear system solve":
+        var rng = initRand(101325)
+        var arr : array[0..(100 * 100 - 1), float64]
+        for i in 0..<arr.len:
+            arr[i] = rng.rand(-100f32..100f32)
+        var mtx  = newHMatrix[float64](100, 100, arr)
+        var lu = mtx.clone()
+        let pivot = lu.lup()
+        var b = newHVector[float64](100, 0.0f)
+        for i in 0...len(b):
+            b[i] = float64(rng.rand(b.len))
+
+        let x = lu.lu_solve(b, pivot)
+        let error  = (mtx * x) - b
+        check(error.norm() < 1e-5)
+    
+    test "triu":
+        let mtx  = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
+        block:
+            let upper = mtx.triu()
+            check(upper == newHMatrix[int8](3, 3, [1i8,2,3,0,5,6,0,0,9]))
+        block:
+            let upper = mtx.triu(1)
+            check(upper == newHMatrix[int8](3, 3, [1i8,2,3,0,1,6,0,0,1]))
+
+    test "tril":
+        let mtx  = newHMatrix[uint](3, 3, [1u,2,3,4,5,6,7,8,9])
+        block:
+            let lower = mtx.tril()
+            check(lower == newHMatrix[uint](3, 3, [1u,0,0,4,5,0,7,8,9]))
+        block:
+            let lower = mtx.tril(1)
+            check(lower == newHMatrix[uint](3, 3, [1u,0,0,4,1,0,7,8,1]))
+
+    test "transpose":
+        block:
+            let mtx = newHMatrix[int](3, 3, [1, 2, 3, 4, 5, 6, 8, 7, 9])
+            let xpose = newHMatrix[int](3, 3, [1,4,8,2,5,7,3,6,9])
+            check(mtx.transpose() == xpose)
+        block:
+            let mtx = newHMatrix[int](2, 5, [1,4,8,2,5,7,3,6,9,0])
+            let xpose = newHMatrix[int](5, 2, [1,7,4,3,8,6,2,9,5,0])
+            check(mtx.transpose() == xpose)

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

@@ -0,0 +1,119 @@
+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, squared_norm2
+from mmath/svector import buildSVector, `-`, abs, norm, Svector
+from mmath/error import SingularMatrixError
+from mmath/rational import newRational, Rational, one, zero, abs, `>`, `<`, `==`, `*`, `+`, `-`, `/`, `/=`, `+=`, `sqrt`
+import unittest
+
+suite "Nim linear algebra library":
+    let test_matrix = newSMatrixFromArray[3, 3, float32]([1.0f32, 2f32, 3f32, 4f32, 5f32, 6f32, 8f32, 7f32, 9f32])
+
+    test "+":
+        let mtx1  = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
+        let mtx2  = newSMatrixFromArray[3, 3, int8]([-1i8,-2,-3,-4,-5,-6,-7,-8,-9])
+        check(mtx1 + mtx2 == SMatrix[3, 3, int8]())
+
+    test "-":
+        let mtx  = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
+        check(mtx - mtx == SMatrix[3, 3, int8]())
+
+    test "+=":
+        var mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
+        mtx += mtx
+        check(mtx == newSMatrixFromArray[3, 3, int8]([2i8,4,6,8,10,12,14,16,18]))
+
+    test "-=":
+        var mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
+        mtx -= mtx
+        check(mtx == SMatrix[3, 3, int8]())
+
+    test "*":
+        block:
+            let mtx  = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,8,7,9])
+            let res = mtx * buildSVector[3, int8](1i8, 2, 3)
+            check(res == buildSVector[3, int8](14i8, 32, 49))
+        block:
+            let mtx = newSMatrixFromArray[2, 5, int]([1,4,8,2,5,7,3,6,9,0])
+            let res = mtx * mtx.transpose()
+            check(res == newSMatrixFromArray[2, 2, int]([110,85,85,175]))
+        block:
+            let mtx = newSMatrixFromArray[2, 5, int]([1,4,8,2,5,7,3,6,9,0])
+            let res = mtx.transpose() * mtx
+            check(res == newSMatrixFromArray[5, 5, int]([
+                50, 25, 50, 65, 5,
+                25, 25, 50, 35, 20,
+                50, 50, 100, 70, 40,
+                65, 35, 70, 85, 10,
+                5, 20, 40, 10, 25]))
+    
+        
+    test "Determinant":
+        check(test_matrix.det() == -9f32)
+
+    test "LU Determinant":
+        check(test_matrix.lu_det() == -9f32)
+
+    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), Rational[int64]]
+        for i in 0..<arr.len:
+            arr[i] = newRational(rng.rand(-20..20).int64, 20.int64)
+        let mtx = newSMatrixFromArray[5, 5, Rational[int64]](arr)
+        var lu = mtx.clone()
+        let pivot = lu.lup()
+        let l = lu.tril(Rational[int64].one)
+        let u = lu.triu()
+        let err = pivot * mtx - (l * u)
+        check(err.squared_norm2() == Rational[int64].zero)
+
+    test "Linear system solve":
+        var rng = initRand(101325)
+        var arr : array[0..(100 * 100 - 1), float64]
+        for i in 0..<arr.len:
+            arr[i] = rng.rand(-100f32..100f32)
+        var mtx  = newSMatrixFromArray[100, 100, float64](arr)
+        var lu = mtx.clone()
+        let pivot = lu.lup()
+        var b : SVector[100, float64]
+        for i in 0...len(b):
+            b[i] = float64(rng.rand(b.len))
+
+        let x = lu.lu_solve(b, pivot)
+        let error  = (mtx * x) - b
+        check(error.norm() < 1e-5)
+    
+    test "triu":
+        let mtx  = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
+        block:
+            let upper = mtx.triu()
+            check(upper == newSMatrixFromArray[3, 3, int8]([1i8,2,3,0,5,6,0,0,9]))
+        block:
+            let upper = mtx.triu(1)
+            check(upper == newSMatrixFromArray[3, 3, int8]([1i8,2,3,0,1,6,0,0,1]))
+
+    test "tril":
+        let mtx  = newSMatrixFromArray[3, 3, uint]([1u,2,3,4,5,6,7,8,9])
+        block:
+            let lower = mtx.tril()
+            check(lower == newSMatrixFromArray[3, 3, uint]([1u,0,0,4,5,0,7,8,9]))
+        block:
+            let lower = mtx.tril(1)
+            check(lower == newSMatrixFromArray[3, 3, uint]([1u,0,0,4,1,0,7,8,1]))
+
+    test "transpose":
+        block:
+            let mtx = newSMatrixFromArray[3, 3, int]([1, 2, 3, 4, 5, 6, 8, 7, 9])
+            let xpose = newSMatrixFromArray[3, 3, int]([1,4,8,2,5,7,3,6,9])
+            check(mtx.transpose() == xpose)
+        block:
+            let mtx = newSMatrixFromArray[2, 5, int]([1,4,8,2,5,7,3,6,9,0])
+            let xpose = newSMatrixFromArray[5, 2, int]([1,7,4,3,8,6,2,9,5,0])
+            check(mtx.transpose() == xpose)
+

vector.nim

@@ -1,72 +0,0 @@
-from oomacro import class
-from utils import `...`
-from sequtils import newSeqWith
-from math import sqrt
-
-type Vector*[T] = seq[T]
-
-#SVector*[S : static[int], T] = array[1..S,T]
-
-proc newVector*[T](size : int, init: T=0) : Vector[T] = 
-    result = newSeq[T](size)
-    for i in 0...len(result):
-        result[i] = init
-
-proc `+`*[T](v1 : Vector[T], v2:Vector[T]) : Vector[T] =
-    result = newVector[T](len(v1))
-    for i in 0...len(v1):
-        result[i] = v1[i] + v2[i]
-
-proc `-`*[T](v1 : Vector[T], v2:Vector[T]) : Vector[T] =
-    result = newVector[T](len(v1))
-    for i in 0...len(v1):
-        result[i] = v1[i] - v2[i]
-
-proc `*`*[T](v1 : Vector[T], v2:Vector[T]) : T =
-    result = 0
-    for i in 0...len(v1):
-        result += v1[i] * v2[i]
-
-proc `+=`*[T](v1 : var Vector[T], v2 : Vector[T]) = 
-    for i in 0...len(v1):
-        v1[i] += v2[i]
-
-proc `-=`*[T](v1 : var Vector[T], v2 : Vector[T]) = 
-    for i in 0...len(v1):
-        v1[i] -= v2[i]
-
-proc `+=`*[T](v : var Vector[T], value : T) = v.add(value)
-
-proc createVector*[T](elems : varargs[T]) : Vector[T] = 
-    result = newSeq[T]()
-    for elem in items(elems):
-        result += elem
-
-proc norm*[T](v : Vector[T]) : T =
-    for value in v:
-        result += v * v
-
-proc abs*[T](v : Vector[T]) : T =
-    return math.sqrt(v.norm)
-
-# let vec = createVector(10,1)
-
-# echo $vec
-# var vec2 = newVector[int](10, 1)
-
-# echo vec2
-
-# type Foo[T, S] = object of RootObj
-#     data : S
-
-# proc `+`[T,S](v1 : Foo[T,S], v2 : Foo[T,S]) : Foo[T,S] =
-#     result = Foo[T,S]()
-#     for i in range(len(v1.data)):
-#         result[i] = v1.data[i] + v2.data[i]
-
-
-# type Bar[T] = Foo[T, array[3,T]]
-
-# var bar = Bar[int]()
-# bar.data[1] = 4
-# bar.data[1] = 4