Algorithm Implementation/Linear Algebra/Determinant of a Matrix

Haskell

--Takes input with the matrix as a list of rows, throws an error if the list isn't square
--This probably isn't the most elegant/efficient way to do it, but it's something
alternatingSum :: Num a => [a] -> a
alternatingSum [] = 0
alternatingSum l  = iter l 0 0
    where iter (x:xs) n result = if even n then iter xs (n + 1) (result + x) else iter xs (n + 1) (result - x)
          iter [] _ result = result

removeColumn :: Int -> [[a]] -> [[a]]
removeColumn n = map (removeIndex n)

removeIndex :: Int -> [a] -> [a]
removeIndex 0 (x:xs) = xs
removeIndex n (x:xs) = x : removeIndex (n-1) xs
removeIndex _ []     = []

isSquare :: [[a]] -> Bool
isSquare (row:rows) = and $ map eqFirstLength rows
    where l = length row
          eqFirstLength list = length list == l
isSquare [] = True

determinant :: Num a => [[a]] -> a
determinant matrix 
    | isSquare matrix   = helper matrix
    | otherwise         = error "DETERMINANT --> MATRIX IS NOT SQUARE"
    where helper [[x1,x2],[y1,y2]] = x1 * y2 - x2 * y1
          helper (x:xs)            = let l = length x in alternatingSum $ zipWith (*) x $ map (\n -> determinant $ removeColumn n xs) [0..l]
          helper [] = 0

Java

public static double[][] reduce(double[][] x, double[][] y, int r, int c, int n) {
    for (int h = 0, j = 0; h < n; ++h) {
        if (h == r) 
            continue;
        for (int i = 0, k = 0; i < n; ++i) {
            if (i == c) 
                continue;
            y[j][k] = x[h][i];
            ++k;
        } //end inner loop
    ++j;
    } //end outer loop
    return y;
} //end method

//===================================================
public static double det(int NMAX, double[][] x) {
    double ret = 0;
    if (NMAX < 4) {//base case
        double prod1 = 1, prod2 = 1;
        for (int i = 0; i < NMAX; ++i) {
            prod1 = 1;
            prod2 = 1;
            for (int j = 0; j < NMAX; ++j) {
                prod1 *= x[(j + i + 1) % NMAX][j];
                prod2 *= x[(j + i + 1) % NMAX][NMAX - j - 1];
            } //end inner loop
            ret += prod1 - prod2;
        } //end outer loop
        return ret;
    } //end base case
    double[][] y = new double [NMAX - 1][NMAX - 1];
    for (int h = 0; h < NMAX; ++h) {
        if (x[h][0] == 0) 
            continue;
        reduce(x, y, h, 0, NMAX);
        if (h % 2 == 0) ret -= det(NMAX - 1, y) * x[h][0];
        if (h % 2 == 1) ret += det(NMAX - 1, y) * x[h][0];
    } //end loop
    return ret;
} //end method

Java

public int determinant(int[][] m) {
	int n = m.length;
	if(n == 1) {
		return m[0][0];
	} else {
		int det = 0;
		for(int j = 0; j < n; j++) {
			det += Math.pow(-1, j) * m[0][j] * determinant(minor(m, 0, j));
		}
		return det;
	}
}

public int[][] minor(final int[][] m, final int i, final int j) {
	int n = m.length;
	int[][] minor = new int[n - 1][n - 1];
	int r = 0, s = 0;
	for(int k = 0; k < n; k++) {
		int[] row = m[k];
		if(k != i) {
			for(int l = 0; l < row.length; l++) {
				if(l != j) {
					minor[r][s++] = row[l];
				}
			}
			r++;
			s = 0;
		}
	}
	return minor;
}

Phix

function minor(sequence a, integer x, integer y)
integer l = length(a)-1
sequence result = repeat(repeat(0,l),l)
    for i=1 to l do
        for j=1 to l do
            result[i][j] = a[i+(i>=x)][j+(j>=y)]
        end for
    end for
    return result
end function
 
function det(sequence a)
    if length(a)=1 then
        return a[1][1]
    end if
    integer sgn = 1
    integer res = 0
    for i=1 to length(a) do
        res += sgn*a[1][i]*det(minor(a,1,i))
        sgn *= -1
    end for
    return res
end function