bool LUDecompose(int n, float *m,
	unsigned short *index, float *detSign)
{
	float *rowNormalizer = new float[n];
	float exchangeParity = 1.0F;
	bool result = false;
	
	// Calculate a normalizer for each row
	for (int i = 0; i < n; i++)
	{
		const float *entry = m + i;
		float maxvalue = 0.0F;
		
		for (int j = 0; j < n; j++)
		{
			float value = fabs(*entry);
			if (value > maxvalue) maxvalue = value;
			entry += n;
		}
		
		if (maxvalue == 0.0F) goto exit; // Singular
		rowNormalizer[i] = 1.0F / maxvalue;
		index[i] = i;
	}
	
	// Perform decomposition
	for (int j = 0; j < n; j++)
	{
		for (int i = 1; i < j; i++)
		{
			// Evaluate Equation (14.15)
			float sum = m[j * n + i];
			for (int k = 0; k < i; k++)
				sum -= m[k * n + i] * m[j * n + k];
			m[j * n + i] = sum;
		}
		
		// Find pivot element
		int pivotRow = -1;
		float maxvalue = 0.0F;
		for (int i = j; i < n; i++)
		{
			// Evaluate Equation (14.17)
			float sum = m[j * n + i];
			for (int k = 0; k < j; k++)
				sum -= m[k * n + i] * m[j * n + k];
			m[j * n + i] = sum;
			
			sum = fabs(sum) * rowNormalizer[i];
			if (sum > maxvalue)
			{
				maxvalue = sum;
				pivotRow = i;
			}
		}
		
		if (pivotRow != j)
		{
			if (pivotRow == -1) goto exit; // Singular
			
			// Exchange rows
			for (int k = 0; k < n; k++)
			{
				float temp = m[k * n + j];
				m[k * n + j] = m[k * n + pivotRow];
				m[k * n + pivotRow] = temp;
			}
			
			unsigned short temp = index[j];
			index[j] = index[pivotRow];
			index[pivotRow] = temp;
			
			rowNormalizer[pivotRow] = rowNormalizer[j];
			exchangeParity = -exchangeParity;
		}
		
		// Divide by pivot element
		if (j != n - 1)
		{
			float denom = 1.0F / m[j * n + j];
			for (int i = j + 1; i < n; i++)
				m[j * n + i] *= denom;
		}
	}
	
	if (detSign) *detSign = exchangeParity;
	result = true;

	exit:
	delete[] rowNormalizer;
	return (result);
}