Determinant.java

package mikera.matrixx.algo;

import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.Matrix33;
import mikera.matrixx.decompose.ILUPResult;
import mikera.matrixx.decompose.impl.lu.SimpleLUP;
import mikera.vectorz.util.IntArrays;

/**
 * Public API function class for determinant calculation.
 * 
 * Normally you should use m.determinant() to calculate the determinant of a matrix. These functions
 * are provided to allow access to alternative algorithms.
 * 
 * @author Mike
 *
 */
public class Determinant {
	
	/**
	 * Calculate the determinant of an AMatrix.
	 * 
	 * @param m
	 * @return
	 */
	public static double calculate(AMatrix m) {
		int rc = m.checkSquare();

		if (rc<=4) {
			// much faster for small matrices
			if (rc<=3) return calculateSmallDeterminant(m,rc);

			// benchmarks show naive method is slightly better for
			// size 4 matrices?
			return naiveDeterminant(m.toMatrix(),rc);
		}
		
		// general determinant uses LUP decomposition
		return calculateLUPDeterminant(m);		
	}
	
	/**
	 * Determinant implemented using the LUP decomposition
	 * 
	 * @param m
	 * @return
	 */
	static double calculateLUPDeterminant(AMatrix m) {
		ILUPResult lup=SimpleLUP.decompose(m);
		double det=lup.getL().diagonalProduct()*lup.getU().diagonalProduct()*lup.getP().determinant();
		return det;
	}
	
	/**
	 * Calculates the determinant of a small square matrix via direct computation
	 * @return
	 */
	static double calculateSmallDeterminant(AMatrix m, int rc) {
		if (rc==1) return m.unsafeGet(0,0);
		if (rc==2) return m.unsafeGet(0,0)*m.unsafeGet(1,1)-m.unsafeGet(1,0)*m.unsafeGet(0,1);
		if (rc==3) {
			return new Matrix33(m).determinant();
		}
		throw new UnsupportedOperationException("Small determinant calculation on size "+rc+" not possible");
	}
	
	/**
	 * Calculate determinant using naive method (brute force)
	 */
	static double naiveDeterminant(Matrix m) {
		return naiveDeterminant(m,m.rowCount());
	}
	
	static double naiveDeterminant(AMatrix m, int rc) {
		int[] inds = new int[rc];
		for (int i = 0; i < rc; i++) {
			inds[i] = i;
		}
		return calcDeterminant(m.toMatrix(),inds, 0);
	}


	private static double calcDeterminant(Matrix m, int[] inds, int offset) {
		int rc = m.rowCount();
		if (offset == (rc - 1))
			return m.unsafeGet(offset, inds[offset]);

		// multiple of first submatrix
		double v0=m.unsafeGet(offset, inds[offset]);
		double det = (v0==0)? 0: v0* calcDeterminant(m,inds, offset + 1);
		
		for (int i = 1; i < (rc - offset); i++) {
			IntArrays.swap(inds, offset, offset + i);
			double v=m.unsafeGet(offset, inds[offset]);
			if (v!=0) {
				det -= v * calcDeterminant(m,inds, offset + 1);
			}
			IntArrays.swap(inds, offset, offset + i);
		}
		return det;
	}

}