Linear algebra

Gauss-Jordan elimination

use math.la: id, scalar

function pivoting(A,B,n,j)
   m = abs(A[j,j])
   k = j
   for i in j+1..n-1
      if m<abs(A[i,j])
         m = abs(A[i,j])
         k = i
      end
   end
   A[j],A[k] = A[k],id(A[j])
   B[j],B[k] = B[k],id(B[j])
end

function gauss_jordan(A,B,n)
   for j in 0..n-1
      pivoting(A,B,n,j)
      B[j] = B[j]/A[j,j]
      A[j] = A[j]/A[j,j]
      for i in j+1..n-1
         if A[i,j]!=0
            B[i] = B[i]/A[i,j]-B[j]
            A[i] = A[i]/A[i,j]-A[j]
         end
      end
   end

   for i in 0..n-2
      for j in i+1..n-1
         B[i] = B[i]-A[i,j]*B[j]
         A[i] = A[i]-A[i,j]*A[j]
      end
   end

   return B
end

# A: a square matrix
# B: a matrix or a vector
function solve(A,B)
   n = A.shape[0]
   A = A.copy
   B = B.copy
   return gauss_jordan(A,B,n)
end

function inv(A)
   n = A.shape[0]
   A = A.copy
   E = scalar(n,1,0)
   return gauss_jordan(A,E,n)
end

Determinant of a matrix

use math.la: id

function pivoting_det(A,n,j)
   m = abs(A[j,j])
   k = j
   for i in j+1..n-1
      if m<abs(A[i,j])
         m = abs(A[i,j])
         k = i
      end
   end
   if k==j
      return false
   else
      A[j],A[k] = A[k],id(A[j])
      return true
   end
end

function det(A)
   n = A.shape[0]
   A = A.copy
   y = 1
   for j in 0..n-1
      if pivoting_det(A,n,j)
         y = -y
      end
      for i in j+1..n-1
         if A[i,j]!=0
            y = y/A[j,j]
            A[i] = A[i]*A[j,j]-A[j]*A[i,j]
         end
      end
      y = y*A[j,j]
   end
   return y
end

Rational numbers

use math.rational: rat
use math.la: matrix
use math.la.inversion: det, inv

A = matrix(
  [1,2,3],
  [4,5,6],
  [7,7,9]
).map(rat)

print(det(A))
# -6

print(inv(A))
# matrix(
#   [-1/2, -1/2, 1/2],
#   [-1, 2, -1],
#   [7/6, -7/6, 1/2]
# )

print(inv(A)*A)
# matrix(
#   [1, 0, 0],
#   [0, 1, 0],
#   [0, 0, 1]
# )

Complex rational numbers

use math.cmath: complex
use math.rational: rat
use math.la: matrix
use math.la.inversion: det, inv

i = complex(0,1)

crat = |z| rat(z.re)+rat(z.im)*i

A = matrix(
  [1+2*i, 3+2*i],
  [2+4*i, 4+3*i]
).map(crat)

print(det(A))
# 0-5*i

print(inv(A))
# matrix(
#   [-3/5+4/5*i, 2/5-3/5*i],
#   [4/5-2/5*i, -2/5+1/5*i]
# )

print(inv(A)*A)
# matrix(
#   [1+0*i, 0+0*i],
#   [0+0*i, 1+0*i]
# )

Random numbers

function random_vector(n,rand)
   return vector(*rand.list(n))
end

function random_matrix(m,n,rand)
   return matrix(*rand.chunks(n).list(m))
end

v = random_vector(4,rng(1..6))
A = random_matrix(4,4,rng(1..6))
B = random_matrix(4,4,rng())