// This is a copy of http://cap-lore.com/MatPhys/determ/d.c
// To compute the determinant of a k by k matrix where A[i, j]
// is located at b + i*s + j.
// k > 0
// the matrix is corrupted.
typedef double R;
R det(R * b, int k, int s){
  if(k==1) return *b;
  {int j; for(j=0; j<k; ++j) if (*(b+s*j) != 0.0) {
     if(j) {int n=k; R * q=b+j*s; while(n--) {
        R t = *(q+n); *(q+n) = *(b+n); *(b+n) = -t;}}
     goto zz;}
  return 0;}
zz: {
  R t = 1./ *b;
  {int i; for(i=1; i<k; ++i) {R * q = b+i*s; R a = t * *q;
     int r; for(r=1; r<k; ++r) *(q+r) -= a * *(b+r);}}
  return *b * det(b+s+1, k-1, s);}}