float det3(float *matrix3)
{
	float det = matrix3[0] * (matrix3[4] * matrix3[8] - matrix3[7] *
							  matrix3[5]) -
			matrix3[1] * (matrix3[3] * matrix3[8] - matrix[6] * matrix[5]) +
			matrix3[2] * (matrix3[3] * matrix3[7] - matrix3[6] * matrix3[4]);
	return det;
}

void submat(float *matrix4, float *matrix3, int i, int j)
{
	int ti, tj, idst, jdst;

	for (ti = 0; ti < 4; ti++)
	{
		if (ti < i)
			idst = ti;
		else if (ti > i)
			idst = ti - 1;
		for (tj = 0; tj < 4; tj++)
		{
			if (tj < j)
				jdst = tj;
			else if (tj > j)
				jdst = tj - 1;
			if (ti != i && tj != j)
				matrix3[idst * 3 + jdst] = matrix4[ti * 4 + tj];
		}
	}
}

float det(float *matrix4)
{
	float det, result = 0;

	i = 1;
	float matrix3[9];
	int n;

	for (n = 0; n < 4; n++, i *= -1)
	{
		submat(mr, matrix3, 0, n);
		det = det3(matrix3);
		result += matrix4[n] * det * i;
	}
	return result;
}

void inverse(float *matrix4a, float *matrix4b)
{
	float mdet = det(matrix4a);
	float temp[9];
	int i, j, sign;

	for (i = 0; i < 4; i++)
		for (j = 0; j < 4; j++)
		{
			sign = 1 - ((i + j) % 2) * 2;
			submat(matrix4b, mtemp, i, j)
					matrix4a[i + j * 4] = (det3(mtemp) * sign) / mdet;
		}
}