added bigint library
This commit is contained in:
253
src/mmath/bigint.nim
Normal file
253
src/mmath/bigint.nim
Normal file
@@ -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
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -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()
|
||||
|
||||
|
||||
@@ -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
|
||||
@@ -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())
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user