libraries/linalg.rt (1312198) - rtensor

linalg.rt

# functions useful in linear algebra (vectors, matrices, and such)
# use this by pasting file contents into the code-area
# matrices are assumed to be row-major

# dot product of u and v
dot(u, v) = if[any((==len(u), ==len(v))), \
    0, \
    u[1] * v[1] + dot(tail(u), tail(v))]

# determinant of matrix m
det(m) = if[len(m) == 1, \
    m[1][1], \
    sum[(k = 1), len(m), \
        (-1)^(k + 1) * m[1][k] * det(minor(m, 1, k))]]

# minor of a m, removing row i and column j
minor(m, i, j) = map(except(m, i), (r => except(r, j)))
except(v, i) = if[all((len(v) == 2, i == 1)), \
    .. v[2], \
    if[len(v) == i, \
        trim(v), \
        (except(trim(v), i), v[len(v)])]]

# cross product of 3D vectors u and v
cross(u, v) = map(1 .. 3, (k => \
    (-1)^(k+1) * det(minor(((0, 0, 0), u, v), 1, k))))

# cross product of 4D vectors u, v, and w
cross(u, v, w) = map(1 .. 4, (k => \
    (-1)^(k+1) * det(minor(((0, 0, 0, 0), u, v, w), 1, k))))

# product of matrix M and vector v
mvm(M, v) = map(v, (_x, i => dot(v, M[i])))

# product of matrices a and b
mm(a, b) = map(1 .. len(a), (i => map(1 .. len(b[1]), (j => \
    dot(a[i], col(b, j))))))

# column i of matrix m
col(m, i) = map(m, (r => r[i]))

# transpose of matrix m
trans(m) = map(1 .. len(m[1]), (k => col(m, k)))

# n-by-n identity matrix
idm(n) = map(1 .. n, (k => zero(k - 1) .. (.. 1) .. zero(n - k)))

# matrix m to the (natural number) power n
mpow(m, n) = mm(m, mpow(m, n - 1))
mpow(m, 0) = idm(len(m))

# matrix m, but in row-echelon form
ref(m) = for[(c = 1), c < len(m[1]), \
    (m = map(m, ((r, i) => if[i < c + 1, \
        r, \
        r - m[c] * r[c] / m[c][c]]))) + (c = c + 1), \
    m]

# internal helper function
nz(x) = if[x == 0, 1, x]

# row-echelon matrix m (answers are probably wrong)
rref(m) = for[(c = 2), c < len(m[1]), \
    (m = map(m, ((r, i) => if[c + 1 < i, \
        r, \
        r - m[c] * r[c] / nz(m[c][c]])))) + (c = c + 1), \
    m]

# invert upper triangular matrix m by recursive partitioning
invut(m) = if[all((len(m) == len(m[1]), \
    all(tail(col(m, 1)) == zero(len(m) - 1)))), \
    if[len(m) == 1, \
        .. .. /m[1][1], \
        (ai, di => \
            ((.. ((.. ai) .. -mm(mm(.. .. ai, .. tail(m[1])), di)[1]))) .. \
            map(tail(m), ((_r, i) => (.. 0) .. di[i])))( \
            /m[1][1], \
            invut(map(tail(m), (r => tail(r)))))], \
    0/0]

rtensor

Import v2.3-dev from a demonic ritual