From a219460d207c2fecf097e8320028f9c085dd7b0c Mon Sep 17 00:00:00 2001
From: Walter Oggioni <oggioni.walter@gmail.com>
Date: Sun, 28 May 2023 13:53:28 +0800
Subject: [PATCH] added bigint library

---
 src/mmath/bigint.nim   | 253 +++++++++++++++++++++++++++++++++++++++++
 src/mmath/hmatrix.nim  | 149 ++++++++++++------------
 src/mmath/hvector.nim  |  28 ++---
 src/mmath/pivot.nim    |   5 +-
 src/mmath/rational.nim |  48 ++++++--
 src/mmath/smatrix.nim  | 133 +++++++++++-----------
 src/mmath/svector.nim  |  39 ++++---
 7 files changed, 475 insertions(+), 180 deletions(-)
 create mode 100644 src/mmath/bigint.nim

diff --git a/src/mmath/bigint.nim b/src/mmath/bigint.nim
new file mode 100644
index 0000000..066b116
--- /dev/null
+++ b/src/mmath/bigint.nim
@@ -0,0 +1,253 @@
+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 argp_parse*(argp : ptr[Argp], argc : cint, argv : cstringArray, flags : cint, arg_index : cint, ctx : pointer) : cint
+    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_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: var 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 `$`*(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
diff --git a/src/mmath/hmatrix.nim b/src/mmath/hmatrix.nim
index 1659825..2e6cb81 100644
--- a/src/mmath/hmatrix.nim
+++ b/src/mmath/hmatrix.nim
@@ -1,8 +1,6 @@
-from sequtils import newSeqWith, map
-from nwo/utils import `...`, `-->`, box
+from nwo/utils import `-->`
 from hvector import HVector, newHVector, buildHVector
 from pivot import HPivot, newHPivot, `[]`, `[]=`, len, SizeError, SingularMatrixError
-from options import Option, none, some
 from math import sqrt
 
 type
@@ -22,27 +20,36 @@ 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:
+    for i in 0..<m.rows:
+        for j in 0..<m.columns:
             yield (i, j, m[i, j])
 
-proc newHMatrix*[T](rows, columns : int) : HMatrix[T] = 
+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 : T) : HMatrix[T] = 
-    result = newHMatrix[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 = newHMatrix[T](rows, columns)
+    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] =
-    result = newHMatrix[T](sz,sz,0)
-    for i in 0...sz:
-        result[i,i] = T.one
+    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 = "["
@@ -61,60 +68,60 @@ proc `$`*[T](m : HMatrix[T]) : string =
 
 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:
+    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, 0)
-    for i in 0...m1.rows:
-        for j in 0...m1.columns:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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] =
@@ -125,8 +132,8 @@ proc `-`*[T](m : HMatrix[T]) : HMatrix[T] =
 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:
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
             if m1[i,j] != m2[i,j]:
                 return false
     return true
@@ -139,7 +146,7 @@ proc transpose*[T](m : HMatrix[T]) : HMatrix[T] =
         result[j, i] = v
 
 proc swap_rows[T](m : var HMatrix[T], id1 : int, id2 : int) =
-    for i in 0...m.columns:
+    for i in 0..<m.columns:
         let tmp = m[id1, i]
         m[id1, i] = m[id2, i]
         m[id2, i] = tmp
@@ -164,7 +171,7 @@ proc add_row[T](
         sourceIndex : int, 
         destIndex : int, 
         factor : T) =
-    for i in 0...m.columns:
+    for i in 0..<m.columns:
         m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor
 
 proc add_row[T](
@@ -180,13 +187,13 @@ 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:
+    for i in 0..<m.rows:
         if m[i, i] == 0:
-            for j in (i + 1)...m.columns:
+            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:
+        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)
@@ -194,13 +201,13 @@ proc gauss_jordan_low*[T](
 proc gauss_jordan_low*[T](
         m : var HMatrix[T]) =
     var pivot = newHPivot[T](m.rows)
-    for i in 0...m.rows:
+    for i in 0..<m.rows:
         if m[i, i] == 0:
-            for j in (i + 1)...m.columns:
+            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:
+        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)
@@ -240,7 +247,7 @@ proc det*[T](m : HMatrix[T]) : T =
     var clone = m.clone()
     clone.gauss_jordan_low()
     result = T.one
-    for i in 0...clone.rows:
+    for i in 0..<clone.rows:
         result *= clone[i, i]
 
 proc invert*[T](m : HMatrix[T]) : HMatrix[T] =
@@ -250,15 +257,15 @@ proc invert*[T](m : HMatrix[T]) : HMatrix[T] =
     result = identity[T](m.rows)
     tmp.gauss_jordan_low(result)
     tmp.gauss_jordan_high(result)
-    for i in 0...result.rows:
+    for i in 0..<result.rows:
         let f = tmp[i, i]
-        for j in 0...result.columns:
+        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:
+    for i in 0..<m.rows:
+        for j in i..<m.columns:
             if i == j:
                 result[i, j] = diag_replace
             else:
@@ -266,15 +273,15 @@ proc triu*[T](m : HMatrix[T], diag_replace: T) : HMatrix[T] =
 
 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:
+    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):
+    for i in 0..<m.rows:
+        for j in 0..<(i + 1):
             if i == j:
                 result[i, j] = diag_replacement
             else:
@@ -282,25 +289,25 @@ proc tril*[T](m : HMatrix[T], diag_replacement : T) : HMatrix[T] =
 
 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):
+    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:
+    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:
+    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:
+    for j in (i + 1)..<m.rows:
         if abs(m[j, i]) > max:
             max = abs(m[i, j])
             max_index = j
@@ -309,24 +316,24 @@ proc lu_pivot[T](m : var HMatrix[T], i : int, pivot : var HPivot[T]) =
 
 proc lup*[T](m : var HMatrix[T]) : HPivot[T] =
     result = newHPivot[T](m.rows)
-    for i in 0...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:
+    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:
+    for i in 0..<m.rows:
         x[i] = b[pivot[i]]
-        for k in 0...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:
+        for k in (i + 1)..<m.rows:
             x[i] = x[i] - m[i, k] * x[k]
-        if m[i,i] != 0:
+        if m[i,i] != T.zero:
             x[i] /= m[i, i]
         else:
             raise newException(SingularMatrixError, "Matrix is singular")
@@ -342,8 +349,8 @@ 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:
+    for i in 0..<m.rows:
+        for j in 0..<m.rows:
             if pivot[j] == i:
                 result[j, i] = T.one
             else:
@@ -362,7 +369,7 @@ proc lu_det*[T](m : var HMatrix[T]) : T =
         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:
+    for i in 0..<m.rows:
         result *= m[i, i]
     if pivot.permutations mod 2 != 0:
         result *= -T.one
@@ -384,7 +391,7 @@ proc norm2*[T](m : HMatrix[T]): T =
 proc `*`*[T](pivot : HPivot[T], m : HMatrix[T]) : HMatrix[T] =
     result = m.clone()
     var pclone = pivot
-    for i in 0...pclone.len():
+    for i in 0..<pclone.len():
         while i != pclone[i]:
             result.swap_rows(i, pclone[i])
             let tmp = pclone[i]
@@ -393,6 +400,6 @@ 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:
+    for i in 0..<pivot.len:
         result[pivot[i], i] = T.one
 
diff --git a/src/mmath/hvector.nim b/src/mmath/hvector.nim
index a37adf6..d89d5ad 100644
--- a/src/mmath/hvector.nim
+++ b/src/mmath/hvector.nim
@@ -1,35 +1,36 @@
-from nwo/utils import `...`
-from sequtils import newSeqWith
-from math import sqrt
-
 type HVector*[T] = seq[T]
 
-proc newHVector*[T](size : int, init: T=0) : HVector[T] = 
+proc newHVector*[T](size : int, init: T=T.zero) : HVector[T] = 
     result = newSeq[T](size)
-    for i in 0...len(result):
+    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):
+    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):
+    for i in 0..<len(v1):
         result[i] = v1[i] - v2[i]
 
 proc `*`*[T](v1 : HVector[T], v2:HVector[T]) : T =
-    result = 0
-    for i in 0...len(v1):
+    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):
+    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):
+    for i in 0..<len(v1):
         v1[i] -= v2[i]
 
 proc `+=`*[T](v : var HVector[T], value : T) = v.add(value)
@@ -40,9 +41,10 @@ proc buildHVector*[T](elems : varargs[T]) : HVector[T] =
         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 math.sqrt(v.norm)
+    return v.norm().sqrt()
 
diff --git a/src/mmath/pivot.nim b/src/mmath/pivot.nim
index bcd142f..9be55b1 100644
--- a/src/mmath/pivot.nim
+++ b/src/mmath/pivot.nim
@@ -1,5 +1,4 @@
 from hvector import HVector, newHVector
-from nwo/utils import `...`
 
 type
     HPivot*[T] = object
@@ -18,7 +17,7 @@ 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:
+    for i in 0..<size:
         result[i] = i
 
 proc `[]`*[SIZE, T](p : SPivot[SIZE, T], index : int) : int = p.data[index]
@@ -26,5 +25,5 @@ proc `[]=`*[SIZE, T](p : var SPivot[SIZE, T], index : int, value : int) = p.data
 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:
+    for i in 0..<SIZE:
         result[i] = i
\ No newline at end of file
diff --git a/src/mmath/rational.nim b/src/mmath/rational.nim
index 5c37af4..fd4856a 100644
--- a/src/mmath/rational.nim
+++ b/src/mmath/rational.nim
@@ -1,8 +1,4 @@
-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 =
+proc gcd*[T](a : T, b : T) : T =
     var n1 = a
     var n2 = b 
     var tmp : T
@@ -12,14 +8,14 @@ proc gcd*[T : SomeInteger](a : T, b : T) : T =
         n2 = tmp mod n2
     return n1
 
-proc mcm*[T : SomeInteger](n1 : T, n2 : T) : T = n1 * n2 div gcd(n1, n2)
+proc mcm*[T](n1 : T, n2 : T) : T = n1 * n2 div gcd(n1, n2)
 
-proc pow*[T : SomeInteger](n : T, e : T) : T =
+proc pow*[T](n : T, e : T) : T =
     result = T.one
     for i in 0..<e:
         result *= n    
 
-type Rational*[T : SomeInteger] = object
+type Rational*[T] = object
     num* : T
     den* : T
 
@@ -34,9 +30,41 @@ proc simplify*[T](self : var Rational[T]) : void =
 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 : int64](_: type[T]) : T = 1.T
 
-proc one*[T](_: typedesc[Rational[T]]): Rational[T] = newRational(T.one, T.one)
+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 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())
 
diff --git a/src/mmath/smatrix.nim b/src/mmath/smatrix.nim
index d31ce04..4f61759 100644
--- a/src/mmath/smatrix.nim
+++ b/src/mmath/smatrix.nim
@@ -1,9 +1,7 @@
-from sequtils import newSeqWith, map
-from nwo/utils import `...`, `-->`, box
+from nwo/utils import `-->`, box
 from svector import SVector
 from pivot import SPivot, newSPivot, `[]`, `[]=`, SingularMatrixError, len
 from math import sqrt
-import random
 
 type
   SMatrix*[ROWS, COLUMNS: static[int], T] = object
@@ -25,19 +23,24 @@ proc `[]=`*[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS,COLUMNS,T], r,
     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:
+    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:
+    for i in 0..<SIZE:
         result[i,i] = T.one
 
 proc `$`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : string = 
@@ -58,60 +61,60 @@ proc `$`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : strin
 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:
+    for i in 0..<result.rows:
+        for j in 0..<result.columns:
             result[i, j] = T.zero
-            for k in 0...m1.columns:
+            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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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:
+    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] =
@@ -119,8 +122,8 @@ proc `-`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatr
         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:
+    for i in 0..<m1.rows:
+        for j in 0..<m1.columns:
             if m1[i,j] != m2[i,j]:
                 return false
     return true
@@ -141,7 +144,7 @@ proc transpose*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) :
         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:
+    for i in 0..<m.columns:
         let tmp = m[id1, i]
         m[id1, i] = m[id2, i]
         m[id2, i] = tmp
@@ -166,7 +169,7 @@ proc add_row[ROWS, COLUMNS : static[int], T](
         sourceIndex : int, 
         destIndex : int, 
         factor : T) =
-    for i in 0...m.columns:
+    for i in 0..<m.columns:
         m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor
 
 proc add_row[ROWS, COLUMNS : static[int], T](
@@ -182,13 +185,13 @@ 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:
+    for i in 0..<m.rows:
         if m[i, i] == T.zero:
-            for j in (i + 1)...m.columns:
+            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:
+        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)
@@ -196,13 +199,13 @@ proc gauss_jordan_low*[ROWS, COLUMNS : static[int], T](
 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:
+    for i in 0..<m.rows:
         if m[i, i] == 0:
-            for j in (i + 1)...m.columns:
+            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:
+        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)
@@ -240,7 +243,7 @@ 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:
+    for i in 0..<clone.rows:
         result *= clone[i, i]
 
 proc invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] =
@@ -248,14 +251,14 @@ proc invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[
     result = identity[SIZE, T]()
     tmp.gauss_jordan_low(result)
     tmp.gauss_jordan_high(result)
-    for i in 0...result.rows:
+    for i in 0..<result.rows:
         let f = tmp[i, i]
-        for j in 0...result.columns:
+        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:
+    for i in 0..<m.rows:
+        for j in 0..<m.columns:
             if i < j:
                 result[i, j] = m[i, j]
             elif i == j:
@@ -264,16 +267,16 @@ proc triu*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replacement:
                 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:
+    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:
+    for i in 0..<m.rows:
+        for j in 0..<m.columns:
             if i > j:
                 result[i, j] = m[i, j]
             elif i == j:
@@ -282,8 +285,8 @@ proc tril*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replacement :
                 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:
+    for i in 0..<m.rows:
+        for j in 0..<m.columns:
             if i >= j:
                 result[i, j] = m[i, j]
             else:
@@ -292,18 +295,18 @@ proc tril*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SI
 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:
+    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:
+    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:
+    for j in (i + 1)..<m.rows:
         if abs(m[j, i]) > max:
             max = abs(m[i, j])
             max_index = j
@@ -311,23 +314,23 @@ proc lu_pivot[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T], i : int, pi
         m.swap_rows(i, max_index, pivot)
 
 proc lu*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) =
-    for i in 0...m.rows:
+    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:
+    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:
+    for i in 0..<m.rows:
         x[i] = b[pivot[i]]
-        for k in 0...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:
+        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]
@@ -342,8 +345,8 @@ proc lu_solve*[SIZE : static[int], 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:
+    for i in 0..<m.rows:
+        for j in 0..<m.rows:
             if pivot[j] == i:
                 result[j, i] = T.one
             else:
@@ -360,7 +363,7 @@ 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 = T.one
-    for i in 0...m.rows:
+    for i in 0..<m.rows:
         result *= m[i, i]
     if pivot.permutations mod 2 != 0:
         result *= -T.one
@@ -372,7 +375,7 @@ proc lu_det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T =
 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():
+    for i in 0..<pclone.len():
         while i != pclone[i]:
             result.swap_rows(i, pclone[i])
             let tmp = pclone[i]
@@ -381,6 +384,6 @@ 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:
+    for i in 0..<SIZE:
         result[pivot[i], i] = T.one
 
diff --git a/src/mmath/svector.nim b/src/mmath/svector.nim
index db3d89c..0a06334 100644
--- a/src/mmath/svector.nim
+++ b/src/mmath/svector.nim
@@ -1,73 +1,76 @@
-from nwo/utils import `...`
-from sequtils import newSeqWith
 from math import sqrt
 
 type SVector*[S : static[int], T] = array[0..(S-1), T]
 
-proc newSVector*[SIZE, T](init: T=0) : SVector[SIZE, T] = 
-    for i in 0...len(result):
+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):
+    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):
+    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):
+    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):
+    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):
+    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):
+    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):
+    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):
+    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):
+    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):
+    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 = 0
-    for i in 0...len(v1):
+    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):
+    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):
+    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