Vector

import church.lang.Array;
import church.lang.ArraySupport;
import church.lang.Interval;

class Vector<T>: vector(T...: elements);

createVector(length, generator) = vector(createArray(length, generator));

createMatrix(rowCount, columnCount, generator) = createVector(rowCount, i −> createVector(columnCount, j −> generator(i, j)));

// Accessors

vector(a).elements = a;

vector(a).size = a.length;

vector(a)[i] = get0(a, i);

// Utilities

commonLength(l1, l2) = {
    if l1 != l2 then error "Vector:: incompatible vectors" else l1;
}

sum1(n, f) = Interval.new(0, n).reduce(additive_identity, s −> i −> s + f(i));

multiplicativeIdentityMatrix(int: N) = {
    createMatrix(N, N, i −> j −> (i == j) ? multiplicative_identity : additive_identity)
}

// Utilities for matrix division / inversion

pivot(A, int: row, int: col) = {
    if (row == A.size) {
        error "Division by singular matrix";
    } else if (A[row][col] != additive_identity) {
        row;
    } else {
        pivot(A, row + 1, col);
    }
}

/**
 * The Gauss-Jordan reduction step, reducing column j of the supplied matrix A to a column vector
 * whose off-diagonal entries are '0' and whose diagonal entry is '1'. Here, '0' and '1' respectively
 * stand for the additive_identity and the multiplicative_identity of the implied element type.
 * The implementation is returned in the form of a matrix operator that may be applied to a general input matrix.
 */
gaussJordanStep(A, int: j) =
    pivot = pivot(A, j, j)
    M −> createVector(A.size, i −> {
        pivotRow = M[pivot]
        if (i == j) {
            d = A[pivot][i]
            createVector(pivotRow.size, k −> pivotRow[k] / d)
        } else {
            source = (i == pivot) ? j : i
            sourceRow = M[source]
            factor = A[source][j] / A[pivot][j]
            createVector(sourceRow.size, k −> sourceRow[k] - factor * pivotRow[k])
        }
    })


gaussJordanElimination(N, D, j) = {
    if (j == D.size) {
        N
    } else {
        op = gaussJordanStep(D, j)
        gaussJordanElimination(op(N), op(D), j + 1);
    }
}

/**
 * Bareiss's 'fraction-free' variant of the gaussJordan reduction step.
 */
bareissStep(A, int: j) =
    pivot = pivot(A, j, j)
    d = j == 0 ? multiplicative_identity : A[j - 1][j - 1]
    M −> createVector(A.size, i −> {
        pivotRow = M[pivot]
        if (i == j) {
//            dd = A[pivot][i]
//            createVector(pivotRow.size, k -> pivotRow[k] / dd)
            pivotRow
        } else {
            source = (i == pivot) ? j : i
            sourceRow = M[source]
            sourceFactor = A[source][j]
            pivotFactor = A[pivot][j]
            createVector(sourceRow.size, k −> (pivotFactor * sourceRow[k] - sourceFactor * pivotRow[k]) /! d)
        }
    })


bareissElimination(rtn, N, D, j) = {
    if (j == D.size) {
        rtn(N, D)
    } else {
        op = bareissStep(D, j)
        bareissElimination(rtn, op(N), op(D), j + 1);
    }
}

cofactors(Vector: m) = bareissElimination(N −> D −> N, multiplicativeIdentityMatrix(m.size), m, 0);

determinant(Vector: m) = bareissElimination(N −> D −> D[m.size - 1][m.size - 1], multiplicativeIdentityMatrix(m.size), m, 0);

// Arithmetic identities / implementation

                                       stream << vector(a) = writeEmbracedElementsToStream(stream, a, "<", ">");

                                    vector(a) == vector(b) = (a == b);

                                 (Vector: a) + (Vector: b) = createVector(commonLength(a.size, b.size), i −> a[i] + b[i]);

                                              -(Vector: a) = createVector(a.size, i −> -a[i]);

                          dotProduct(Vector: a, Vector: b) = sum1(commonLength(a.size, b.size), i −> a[i] * b[i]);

// Matrix * Vector
//                                 (Vector: m) * (Vector: v) = createVector(m.size, i -> dotProduct(m[i], v));

// Matrix * Matrix
                                 (Vector: m) * (Vector: n) = createVector(m.size, i −> createVector(n[i].size,
                                                                                      j −> sum1(m[i].size,
                                                                                          k −> m[i][k] * n[k][j])));

// Matrix / Matrix
                                 (Vector: n) / (Vector: d) = gaussJordanElimination(n, d, 0);

// Matrix inversion
                                              /(Vector: d) = gaussJordanElimination(
                                                                multiplicativeIdentityMatrix(d.size), d, 0);