summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lalge/README.md46
-rw-r--r--lalge/src/matrices.nim171
-rw-r--r--lalge/src/types.nim8
-rw-r--r--lalge/src/vectors.nim57
-rw-r--r--lalge/tests/nim.cfg1
-rw-r--r--lalge/tests/tests.nim144
6 files changed, 427 insertions, 0 deletions
diff --git a/lalge/README.md b/lalge/README.md
new file mode 100644
index 0000000..640654d
--- /dev/null
+++ b/lalge/README.md
@@ -0,0 +1,46 @@
+# lalge (name tbd)
+
+A functional, dead simple linear algebra library written in Nim. Mostly for learning purposes.
+
+This is also serving as my finals review, haha.
+
+If you're looking for an actually robust linear algebra library for Nim, check out [Neo](https://github.com/andreaferretti/neo).
+
+## What's this language?
+
+Nim is a compiled whitespace-sensitive systems language with a strong type system and interesting memory management capabilities.
+
+It doesn't quite beat Python's usability or Rust's features, but I like it quite a lot, and think it's a good complement to other languages I use.
+
+Nim's [UFCS](https://en.wikipedia.org/wiki/Uniform_Function_Call_Syntax) (TL;DR: `a.add(b) == add(a, b)`) and operator overloading support is also a major draw for this type of project.
+
+## Architecture
+
+lalge uses Nim's `seqs` - dynamic arrays, equivalent to Rust's `Vec`, Go's `slice`, or Java's `ArrayList` - to implement matrices. This may have some performance overhead. On the other hand, it has the benefit of making the code much simpler.
+
+This also has the benefit of being able to declare matrices as such:
+```nim
+example: Matrix = @[
+ @[1, 0, 0],
+ @[0, 1, 0],
+ @[0, 0, 1],
+]
+```
+
+## Goals
+- [ ] Cover the basic operations associated with the MATH 221 / typical introductory linear algebra curriculum
+ - [x] Addition / Subtraction
+ - [x] (Scalar / Vector) Multiplication
+ - [x] Inverse
+ - [x] Transpose
+ - [ ] Eigenvalues / Eigenvectors
+ - [ ] everything that i am forgetting
+- [ ] Helpful characteristic functions
+ - [ ] Dimension
+ - [ ] Rank / Nullity
+ - [ ] Spans
+ - [ ] Determinant
+- [x] Support arbitrary-sized matrices with internal seq representation
+- [x] Support specialized "vectors", as opposed to Nx1 matrices
+
+After I get to a reasonable point of completeness, I'd love to go through and compare my library to Eigen (generally praised as a Very Good library)
diff --git a/lalge/src/matrices.nim b/lalge/src/matrices.nim
new file mode 100644
index 0000000..8e3e680
--- /dev/null
+++ b/lalge/src/matrices.nim
@@ -0,0 +1,171 @@
+import std/math, std/sugar
+import types
+
+# A few notes about Nim:
+
+# "result" lets us operate on the return type
+# without (necessarily) constructing it.
+# this is very useful, and as such i use it a lot.
+
+# nim also supplies many functional tools which
+# i also made liberal use of, occasionally to the
+# detriment of readability. sorry.
+
+# the "*" operator on procs is the public marker.
+# oh, and functions with "side effects" are called "procs".
+
+# TODO: implement a better error handling system
+const not_square = "The matrix is not square"
+const mismatch = "Mismatched rows and columns"
+
+## Overrides the echo() proc to provide more human-readable output
+proc echo*(A: Matrix) =
+ echo "@["
+ for row in A:
+ echo " ", row
+ echo "]"
+
+## Generate a new, filled XxX Matrix
+func gen*(rows, columns: Natural, element: float = 0.0): Matrix =
+ result.setLen(rows)
+ for i in 0 ..< rows:
+ result[i].setLen(columns)
+ for j in 0 ..< columns:
+ result[i][j] = element
+
+## Returns the number of rows in a matrix
+func rows*(A: Matrix): int =
+ return A.len()
+
+## Returns the number of columns in a matrix
+func cols*(A: Matrix): int =
+ if A.rows() == 0:
+ return 0
+ else:
+ return A[0].len()
+
+## Simple rows iterator for uniformity
+iterator rows*(A: Matrix): RowVector =
+ for row in A:
+ yield row
+
+## Simple column iterator for matrices
+iterator cols*(A: Matrix): Vector =
+ for i in 0 ..< A.cols():
+ var row: Vector
+ for j in 0 ..< A.rows():
+ row.add(A[i][j])
+ yield row
+
+## Apply an arbitrary function to every element of a matrix
+func map*(A: Matrix, oper: proc(i, j: int): float): Matrix =
+ result = gen(A.rows(), A.cols())
+ for i in 0 ..< A.rows():
+ for j in 0 ..< A.cols():
+ # Note that the position is _actively_ provided
+ result[i][j] = oper(i, j)
+
+# TODO: check if these meet assert requirements
+
+## Matrix addition
+func `+`*(A, B: Matrix): Matrix =
+ return map(A, (i, j) => A[i][j] + B[i][j])
+
+## Matrix subtraction
+func `-`*(A, B: Matrix): Matrix =
+ return map(A, (i, j) => A[i][j] - B[i][j])
+
+## Scalar-Matrix multiplication
+func `*`*(a: float, B: Matrix): Matrix =
+ return map(B, (i, j) => a * B[i][j])
+
+## Scalar-Matrix multiplication
+# func `*`*(A: Matrix, b: int): Matrix =
+# return b * A
+
+## Matrix multiplication
+func `*`*(A, B: Matrix): Matrix =
+ assert A.rows() == B.cols(), mismatch
+ result = gen(A.rows(), B.cols())
+ for i in 0 ..< A.rows():
+ for j in 0 ..< B.cols():
+ for k in 0 ..< B.rows():
+ result[i][j] += A[i][k] * B[k][j]
+
+## Absolute value of a matrix
+func abs*(A: Matrix): Matrix =
+ return map(A, (i, j) => abs(A[i][j]))
+
+## Simple row function for uniformity
+func row*(A: Matrix, r: int): RowVector =
+ return A[r]
+
+## Returns a "column vector" of a matrix as a Xx1 Matrix
+# func col*(A: Matrix, j: int): Vector =
+# result = gen(A.rows(), 1)
+# for i in 0 ..< A.rows():
+# result[i] = @[A[i][j]]
+
+## Alternate implementation of col, returns a Vector
+func col*(A: Matrix, c: int): Vector =
+ for row in A:
+ result.add(@[row[c]])
+
+## Produce a vector of the diagonal entries of a matrix
+func diag*(A: Matrix): Vector =
+ assert A.rows() == A.cols(), not_square
+ for i in 0 ..< A.rows():
+ result.add(A[i][i])
+
+## Transpose the rows and columns of a matrix
+func transpose*(A: Matrix): Matrix =
+ result = gen(A.cols(), A.rows())
+ return map(result, (i, j) => (A[j][i]))
+
+# func transpose*(a: Vector): Matrix =
+# result = gen(a.len(), 1)
+# return map(result, (i, j) => (a[j]))
+
+## Generate an arbitary sized identity matrix
+func identity*(size: Natural): Matrix =
+ return map(gen(size, size), (i, j) => (if i==j: 1.0 else: 0.0))
+
+## Oft-used identity matrices
+let
+ I1*: Matrix = identity(1)
+ I2*: Matrix = identity(2)
+ I3*: Matrix = identity(3)
+ I4*: Matrix = identity(4)
+ I5*: Matrix = identity(5)
+
+## Calculates the determinant of a matrix through Laplace expansion
+func det*(A: Matrix): float =
+ assert A.rows() == A.cols(), not_square
+ # Shortcut by formula
+ if A.rows() == 2:
+ return (A[0][0]*A[1][1]) - (A[0][1]*A[1][0])
+ # Actual factual formula
+ for i in 0 ..< A.rows():
+ var sub: Matrix
+ for a in 1 ..< A.rows():
+ var row: RowVector
+ for b in 0 ..< A.cols():
+ if b == i: continue
+ row.add(A[a][b])
+ sub.add(row)
+ result += (-1.0)^i * A[0][i] * det(sub)
+
+# func spans(A, B: Matrix): bool =
+# return false
+# func subspace(a, b: Matrix): bool =
+# return false
+
+# Matrix Properties
+
+## Calculate the dimension of a matrix by Gaussian Elimination
+# func dim(a: Matrix): int =
+# return 1
+# func rank(a: Matrix): int =
+# return 1
+# func nullity(a: Matrix): int =
+# return 1
diff --git a/lalge/src/types.nim b/lalge/src/types.nim
new file mode 100644
index 0000000..321042a
--- /dev/null
+++ b/lalge/src/types.nim
@@ -0,0 +1,8 @@
+type
+ # Any generic Vector is assumed to be a column vector.
+ Vector* = seq[float]
+ RowVector* = Vector # can this be distinct?
+ Matrix* = seq[RowVector]
+
+type
+ Radian* = float
diff --git a/lalge/src/vectors.nim b/lalge/src/vectors.nim
new file mode 100644
index 0000000..a9880d9
--- /dev/null
+++ b/lalge/src/vectors.nim
@@ -0,0 +1,57 @@
+import std/math, std/sequtils, std/sugar
+import types
+
+const dim_mismatch = "Dimensional mismatch - check your vectors"
+
+## Vector addition
+func `+`*(a, b: Vector): Vector =
+ assert a.len() == b.len(), dim_mismatch
+ result.setLen(a.len)
+ for i in 0 ..< a.len():
+ result[i] = (a[i] + b[i])
+
+## Vector subtraction
+func `-`*(a, b: Vector): Vector =
+ assert a.len() == b.len(), dim_mismatch
+ result.setLen(a.len)
+ for i in 0 ..< a.len():
+ result[i] = (a[i] - b[i])
+
+## Vector dot product
+func `*`*(a, b: Vector): float =
+ assert a.len() == b.len(), dim_mismatch
+ for i in 0 ..< a.len():
+ result += a[i] * b[i]
+
+func dot*(a, b: Vector): float =
+ return a * b
+
+## Scalar-Vector multiplication
+func `*`*(a: float, b: Vector): Vector =
+ return map(b, (x) => (a*x))
+
+## Produce the length (in space) of a vector
+func length*(a: Vector): float =
+ return sqrt(a * a)
+
+## Produce the number of elements in a vector
+func size*(a: Vector): int =
+ return len(a)
+
+## Returns whether two vectors are orthagonal
+func ortho*(a, b: Vector): bool =
+ return a * b == 0
+
+## Produce the angle between two vectors, in radians
+func angle*(a, b: Vector): Radian =
+ return arccos((a * b) / (a.length * b.length))
+
+## Produce the cross product between two 3D vectors
+func cross*(a, b: Vector): Vector =
+ assert a.len() == 3, "The first vector is not three-dimensional"
+ assert b.len() == 3, "The second vector is not three-dimensional"
+ return @[
+ (a[1]*b[2]) - (b[1]*a[2]),
+ (a[2]*b[0]) - (b[2]*a[0]),
+ (a[0]*b[1]) - (b[0]*a[1])
+ ]
diff --git a/lalge/tests/nim.cfg b/lalge/tests/nim.cfg
new file mode 100644
index 0000000..85bf6c4
--- /dev/null
+++ b/lalge/tests/nim.cfg
@@ -0,0 +1 @@
+--path:"../src/"
diff --git a/lalge/tests/tests.nim b/lalge/tests/tests.nim
new file mode 100644
index 0000000..e798c39
--- /dev/null
+++ b/lalge/tests/tests.nim
@@ -0,0 +1,144 @@
+import std/math
+import types, vectors, matrices
+
+let A: Matrix = @[
+ @[10.0, 2.0, 3.0, 5.0, 6.0],
+ @[8.0, 7.0, 6.0, 4.0, 3.0],
+ @[4.0, 6.0, 0.0, 9.0, 0.0],
+ @[6.0, 7.0, 9.0, 3.0, 9.0],
+ @[3.0, 0.0, 7.0, 9.0, 9.0],
+]
+let B: Matrix = @[
+ @[1.0, 2.0, 3.0, 4.0, 5.0],
+ @[5.0, 3.0, 1.0, 22.0, 3.0],
+ @[5.0, 21.0, 4.0, 6.0, 3.0],
+ @[12.0, 1.0, 5.0, 0.0, 9.0],
+ @[6.0, 7.0, 1.0, 3.0, 5.0],
+]
+
+let L: Matrix = @[
+ @[1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+ @[-1.0, 1.0, 0.0, 0.0, 0.0, 0.0],
+ @[0.0, -1.0, 1.0, 0.0, 0.0, 0.0],
+ @[0.0, 0.0, -1.0, 1.0, 0.0, 0.0],
+ @[0.0, 0.0, 0.0, -1.0, 1.0, 0.0],
+ @[0.0, 0.0, 0.0, 0.0, -1.0, 1.0],
+]
+let U: Matrix = @[
+ @[1.0, -1.0, 0.0, 0.0, 0.0, 0.0],
+ @[0.0, 1.0, -1.0, 0.0, 0.0, 0.0],
+ @[0.0, 0.0, 1.0, -1.0, 0.0, 0.0],
+ @[0.0, 0.0, 0.0, 1.0, -1.0, 0.0],
+ @[0.0, 0.0, 0.0, 0.0, 1.0, -1.0],
+ @[0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+]
+
+let
+ P13: Matrix = @[
+ @[0.0, 0.0, 1.0],
+ @[0.0, 1.0, 0.0],
+ @[1.0, 0.0, 0.0],
+ ]
+ P23: Matrix = @[
+ @[1.0, 0.0, 0.0],
+ @[0.0, 0.0, 1.0],
+ @[0.0, 1.0, 0.0],
+ ]
+
+let foo = @[
+ @[1.0, 1.0, 1.0],
+ @[3.0, 4.0, 5.0],
+ @[40.0, 51.0, 12.0]
+]
+
+let eliminationMatrix: Matrix = @[
+ @[1.0, 0.0, 0.0],
+ @[-2.0, 1.0, 0.0],
+ @[0.0, 0.0, 1.0]
+]
+
+let
+ ones: Matrix = @[
+ @[1.0, 2.0, 3.0],
+ @[4.0, 5.0, 6.0]
+ ]
+ twos: Matrix = @[
+ @[7.0, 8.0],
+ @[9.0, 10.0],
+ @[11.0, 12.0]
+ ]
+
+# Addition and subtraction
+assert I2 + I2 == 2 * I2
+assert I3-I3 == gen(3, 3, 0.0)
+
+# Num. of rows / columns
+assert rows(I2) == 2
+assert I2.rows() == 2
+assert I2.rows == 2
+assert cols(I3) == 3
+assert I3.cols() == 3
+assert I3.cols == 3
+
+# Identity Matrices
+assert I1 == @[
+ @[1.0],
+]
+assert I2 == @[
+ @[1.0, 0.0],
+ @[0.0, 1.0],
+]
+assert I3 == @[
+ @[1.0, 0.0, 0.0],
+ @[0.0, 1.0, 0.0],
+ @[0.0, 0.0, 1.0],
+]
+assert I4 == @[
+ @[1.0, 0.0, 0.0, 0.0],
+ @[0.0, 1.0, 0.0, 0.0],
+ @[0.0, 0.0, 1.0, 0.0],
+ @[0.0, 0.0, 0.0, 1.0],
+]
+assert I5 == @[
+ @[1.0, 0.0, 0.0, 0.0, 0.0],
+ @[0.0, 1.0, 0.0, 0.0, 0.0],
+ @[0.0, 0.0, 1.0, 0.0, 0.0],
+ @[0.0, 0.0, 0.0, 1.0, 0.0],
+ @[0.0, 0.0, 0.0, 0.0, 1.0],
+]
+
+assert det(@[
+ @[2.0, -1.0, 0.0, -1.0],
+ @[-1.0, 2.0, -1.0, 0.0],
+ @[0.0, -1.0, 2.0, -1.0],
+ @[-1.0, 0.0, -1.0, 2.0],
+]) == 0
+
+assert det(A * B) == (det(A) * det(B))
+assert det(I3) == 1
+assert det(2*I4) == 16
+
+# Vectors
+
+let
+ alice: Vector = @[1.0, 1.0, 1.0]
+ bob = @[3.0, 4.0, 5.0]
+ carol = @[40.0, 51.0, 12.0]
+ dan: Vector = @[10.0, 20.0, 30.0, 40.0]
+ eve = @[2.0, 3.0, 4.0]
+
+assert alice + bob == bob + alice
+assert alice + carol == carol + alice
+assert carol - eve != eve - carol
+assert alice.length() == sqrt(3.0)
+
+# echo transpose(ones)
+# echo P13 * 8
+# echo P13 + P23
+# echo ones * twos
+# assert col(eliminationMatrix, 0) == col(eliminationMatrix, 0)
+# echo identity(0)
+# echo identity(2)
+# echo (1/0)
+# echo L*U
+# echo U*L