How do I calculate the determinant of a 6x6 matrix in C++ really fast?

Question

In C++ I need to calculate the determinant of a 6x6 matrix really fast.

This is how I would do this for a 2x2 matrix:

double det2(double A[2][2]) {
   return A[0][0]*A[1][1] - A[0][1]*A[1][0];
}

I want a similar function for the determinant of a 6x6 matrix but I do not want to write it by hand since it contains 6! = 720 terms where each term is the product of 6 elements in the matrix.

Therefore I want to use Leibniz formula:

static int perms6[720][6];
static int signs6[720];
double det6(double A[6][6]) {
  double sum = 0.0;
  for(int i = 0; i < 720; i++) {
     int j0 = perms6[i][0];
     int j1 = perms6[i][1];
     int j2 = perms6[i][2];
     int j3 = perms6[i][3];
     int j4 = perms6[i][4];
     int j5 = perms6[i][5];
     sum += signs6[i]*A[0]*A[j0]*A[1]*A[j1]*A[2]*A[j2]*A[3]*A[j3]*A[4]*A[j4]*A[5]*A[j5];
  }
  return sum; 
}

How do I find the permutations and the signs?

Is there some way I could get the compiler to do more of the work (e.g. C macros or template metaprogramming) so that the function would be even faster?

EDIT: I just timed the following code (Eigen):

Matrix<double,6,6>  A;
// ... fill A

for(long i = 0; i < 1e6; i++) {
    PartialPivLU< Matrix<double,6,6> > LU(A);
    double d = LU.determinant();
}

to 1.25 s. So using LU or Gauss decomposition is definitely fast enough for my use!

Answer 1

Use Gauss method to make the matrix upper-triangular. For every operation you know how determinant is changed (not changed of multiplied by constant d) and it works in O(n^3). After that just multiply numbers on main diagonal and delete to product of all d's