/* This is a quick-n-dirty implementation for computing the determinant
 * of a matrix.
 */

#define MS 100

typedef long double matrix[MS][MS];
typedef long double vector[MS];


typedef long double my_dbl;

/* alternate() returns the determinant of a submatrix of matrix m.
 *   vlen: input, size of submatrix within m
 *   ir, ic: vectors containing row numbers and column numbers
 *         from which the alternating product should be composed.
 *
 *  The implementation is recursive: the alternating product for
 *  larger matrices is expressed in terms of smaller ones.  The
 *  cases for the first 5 orders are computed explicitly, to minimize
 *  the overhead of a recursive call for these smaller sizes.  
 *
 *  The code for the first 5 cases was generated automatically
 *  to ensure correctness, i.e. to avoid typing mistakes.  See
 *  the ifdef'd portion below for the generator.
 *  
 *  Note that execution time grows like N factorial. Thus n=14
 *  might take less than an hour, n=15 might take more, and n=16 
 *  might take days.
 */
my_dbl
alternate (matrix *m, vector *ir, vector *ic, int vlen)
{
	my_dbl det;
	if (2 == vlen)
	{
		int ra, rb, ca, cb;
		ra = (*ir)[0];
		rb = (*ir)[1];

		ca = (*ic)[0];
		cb = (*ic)[1];

		det = (*m)[ra][ca] * (*m)[rb][cb];
		det -= (*m)[ra][cb] * (*m)[rb][ca];
		return det;
	}
	else if (3 == vlen)
	{
		int r0, r1, r2, c0, c1, c2;
		r0 = (*ir)[0];
		r1 = (*ir)[1];
		r2 = (*ir)[2];

		c0 = (*ic)[0];
		c1 = (*ic)[1];
		c2 = (*ic)[2];

		det = 
		(  ((*m)[r0][c0]) * (((*m)[r1][c1]) * ((*m)[r2][c2]) - ((*m)[r1][c2]) * ((*m)[r2][c1]))
		 - ((*m)[r0][c1]) * (((*m)[r1][c0]) * ((*m)[r2][c2]) - ((*m)[r1][c2]) * ((*m)[r2][c0]))
		 + ((*m)[r0][c2]) * (((*m)[r1][c0]) * ((*m)[r2][c1]) - ((*m)[r1][c1]) * ((*m)[r2][c0]))
		);
		return det;
	}
	else if (4 == vlen)
	{
		int r0, r1, r2, r3, c0, c1, c2, c3;
		r0 = (*ir)[0];
		r1 = (*ir)[1];
		r2 = (*ir)[2];
		r3 = (*ir)[3];

		c0 = (*ic)[0];
		c1 = (*ic)[1];
		c2 = (*ic)[2];
		c3 = (*ic)[3];

		det = 
		(  ((*m)[r0][c0]) * (  ((*m)[r1][c1]) * (((*m)[r2][c2]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c2]))
		 - ((*m)[r1][c2]) * (((*m)[r2][c1]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c1]))
		 + ((*m)[r1][c3]) * (((*m)[r2][c1]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c1]))
		)
		 - ((*m)[r0][c1]) * (  ((*m)[r1][c0]) * (((*m)[r2][c2]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c2]))
		 - ((*m)[r1][c2]) * (((*m)[r2][c0]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c0]))
		 + ((*m)[r1][c3]) * (((*m)[r2][c0]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c0]))
		)
		 + ((*m)[r0][c2]) * (  ((*m)[r1][c0]) * (((*m)[r2][c1]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c1]))
		 - ((*m)[r1][c1]) * (((*m)[r2][c0]) * ((*m)[r3][c3]) - ((*m)[r2][c3]) * ((*m)[r3][c0]))
		 + ((*m)[r1][c3]) * (((*m)[r2][c0]) * ((*m)[r3][c1]) - ((*m)[r2][c1]) * ((*m)[r3][c0]))
		)
		 - ((*m)[r0][c3]) * (  ((*m)[r1][c0]) * (((*m)[r2][c1]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c1]))
		 - ((*m)[r1][c1]) * (((*m)[r2][c0]) * ((*m)[r3][c2]) - ((*m)[r2][c2]) * ((*m)[r3][c0]))
		 + ((*m)[r1][c2]) * (((*m)[r2][c0]) * ((*m)[r3][c1]) - ((*m)[r2][c1]) * ((*m)[r3][c0]))
		));
		return det;
	}
	else if (5 == vlen)
	{
		int r0, r1, r2, r3, r4, c0, c1, c2, c3, c4;
		r0 = (*ir)[0];
		r1 = (*ir)[1];
		r2 = (*ir)[2];
		r3 = (*ir)[3];
		r4 = (*ir)[4];

		c0 = (*ic)[0];
		c1 = (*ic)[1];
		c2 = (*ic)[2];
		c3 = (*ic)[3];
		c4 = (*ic)[4];

		det = 
		(  ((*m)[r0][c0]) * (  ((*m)[r1][c1]) * (  ((*m)[r2][c2]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
		 - ((*m)[r2][c3]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
		)
		 - ((*m)[r1][c2]) * (  ((*m)[r2][c1]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
		 - ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
		)
		 + ((*m)[r1][c3]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
		)
		 - ((*m)[r1][c4]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
		 + ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
		)
		)
		 - ((*m)[r0][c1]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c2]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
		 - ((*m)[r2][c3]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
		)
		 - ((*m)[r1][c2]) * (  ((*m)[r2][c0]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
		 - ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
		)
		 + ((*m)[r1][c3]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
		)
		 - ((*m)[r1][c4]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
		)
		)
		 + ((*m)[r0][c2]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c1]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
		 - ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
		)
		 - ((*m)[r1][c1]) * (  ((*m)[r2][c0]) * (((*m)[r3][c3]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c3]))
		 - ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
		)
		 + ((*m)[r1][c3]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
		 - ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
		)
		 - ((*m)[r1][c4]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
		 - ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
		)
		)
		 - ((*m)[r0][c3]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
		)
		 - ((*m)[r1][c1]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
		)
		 + ((*m)[r1][c2]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c1]))
		 - ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c4]) - ((*m)[r3][c4]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c4]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
		)
		 - ((*m)[r1][c4]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
		 - ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
		)
		)
		 + ((*m)[r0][c4]) * (  ((*m)[r1][c0]) * (  ((*m)[r2][c1]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
		 + ((*m)[r2][c3]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
		)
		 - ((*m)[r1][c1]) * (  ((*m)[r2][c0]) * (((*m)[r3][c2]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c2]))
		 - ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
		)
		 + ((*m)[r1][c2]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c1]))
		 - ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c3]) - ((*m)[r3][c3]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c3]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
		)
		 - ((*m)[r1][c3]) * (  ((*m)[r2][c0]) * (((*m)[r3][c1]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c1]))
		 - ((*m)[r2][c1]) * (((*m)[r3][c0]) * ((*m)[r4][c2]) - ((*m)[r3][c2]) * ((*m)[r4][c0]))
		 + ((*m)[r2][c2]) * (((*m)[r3][c0]) * ((*m)[r4][c1]) - ((*m)[r3][c1]) * ((*m)[r4][c0]))
		)));
		return det;
	}
	else
	{
		vector row, col;

		int rrr = (*ir)[0];
		
		// skip the top row when working submatrices
		int i,j;
		for (i=0; i<vlen-1; i++)
		{
			row[i] = (*ir)[i+1];
		}
		
		my_dbl sign = 1.0L;
		my_dbl acc = 0.0L;
		for (i=0; i<vlen; i++)
		{
			int k=0;
			for (j=0; j<vlen-1; j++)
			{
				if (k==i) k++;  // skip this column
				col[j] = (*ic)[k];
				k++;
			}

			int ccc = (*ic)[i];
			acc += sign * (*m)[rrr][ccc] * alt (m, &row, &col, vlen-1);
			sign = -sign;
		}
		return acc;
	}
}

/* Return determinant of matrix m, of dimension dim */

my_dbl
determinant (matrix *m, int dim)
{
	int i;
	vector row, col;

	for (i=0; i<dim; i++)
	{
		row[i] = i;
		col[i] = i;
	}

	return alternate (m, &row, &col, dim);
}

#ifdef TEST
main ()
{
	matrix m;
	int  i, j;
	for (i=0; i<5; i++) 
	{
		for (j=0; j<5; j++)
		{
			// m[i][j] = j*i;
			m[i][j] = 0.0;
		}
		m[i][i] = 1.0;
	}

	printf ("its %Lg\n", determinant (&m, 5));
	return 0;
}
#endif

/* ===================================================================== */
#ifdef DET_GENERATOR

typedef int vector[30];

/* This code generates the source code for the lower orders.
 * Its straightforward to audit
 */
void
alt (vector *ir, vector *ic, int vlen)
{
	if (2 == vlen)
	{
		int ra, rb, ca, cb;
		ra = (*ir)[0];
		rb = (*ir)[1];

		ca = (*ic)[0];
		cb = (*ic)[1];

		printf ("(((*m)[r%d][c%d]) * ((*m)[r%d][c%d]) - ((*m)[r%d][c%d]) * ((*m)[r%d][c%d]))\n",
			ra,ca,rb,cb, ra,cb,rb,ca);
	}
	else
	{
		vector row, col;

		int rrr = (*ir)[0];
		int i,j;
		for (i=0; i<vlen-1; i++)
		{
			row[i] = (*ir)[i+1];
		}

		int sign = 1;
		printf ("(  ");
		for (i=0; i<vlen; i++)
		{
			int k=0;
			for (j=0; j<vlen-1; j++)
			{
				if (k==i) k++;  // skip this column
				col[j] = (*ic)[k];
				k++;
			}

			int ccc = (*ic)[i];
			if (i != 0) {
				if (sign > 0) printf (" + ");
				else printf (" - ");
			}
			printf ("((*m)[r%d][c%d]) * ", rrr, ccc);
			alt (&row, &col, vlen-1);
			sign = -sign;
		}
		printf (")\n");
	}
}

void
det_gen (int dim)
{
	int i;
	vector row, col;

	for (i=0; i<dim; i++)
	{
		row[i] = i;
		col[i] = i;
	}

	alt (&row, &col, dim);
}

main ()
{
	det_gen (5);
	return 0;
}
#endif