Calculating an NxN matrix determinant in C#

How do you calculate the determinant of an NxN ma开发者_JAVA百科trix C# ?

The OP posted another question asking specifically about 4x4 matrices, which has been closed as an exact duplicate of this question. Well, if you're not looking for a general solution but instead are constrained to 4x4 matrices alone, then you can use this ugly looking but tried-and-true code:

public double GetDeterminant() {
    var m = _values;
    return
         m[12] * m[9]  * m[6]  * m[3]   -  m[8] * m[13] * m[6]  * m[3]   -
         m[12] * m[5]  * m[10] * m[3]   +  m[4] * m[13] * m[10] * m[3]   +
         m[8]  * m[5]  * m[14] * m[3]   -  m[4] * m[9]  * m[14] * m[3]   -
         m[12] * m[9]  * m[2]  * m[7]   +  m[8] * m[13] * m[2]  * m[7]   +
         m[12] * m[1]  * m[10] * m[7]   -  m[0] * m[13] * m[10] * m[7]   -
         m[8]  * m[1]  * m[14] * m[7]   +  m[0] * m[9]  * m[14] * m[7]   +
         m[12] * m[5]  * m[2]  * m[11]  -  m[4] * m[13] * m[2]  * m[11]  -
         m[12] * m[1]  * m[6]  * m[11]  +  m[0] * m[13] * m[6]  * m[11]  +
         m[4]  * m[1]  * m[14] * m[11]  -  m[0] * m[5]  * m[14] * m[11]  -
         m[8]  * m[5]  * m[2]  * m[15]  +  m[4] * m[9]  * m[2]  * m[15]  +
         m[8]  * m[1]  * m[6]  * m[15]  -  m[0] * m[9]  * m[6]  * m[15]  -
         m[4]  * m[1]  * m[10] * m[15]  +  m[0] * m[5]  * m[10] * m[15];
}

It assumes you store your vector data in a 16-element array called _values (of double in this case, but float would work too), in the following order:

0, 1, 2, 3,
4, 5, 6, 7,
8, 9, 10, 11,
12, 13, 14, 15

Reduce to upper triangular form, then make a nested loop where you multiply all the values at position i == j together. There you have it.

The standard method is LU decomposition. You may want to use a library instead of coding it yourself. I don't know about C#, but the 40-year standard is LAPACK.

This solution is achieved using row operations. I took an identity matrix of the same dimensions as of my target matrix and then converted the target matrix to the identity matrix such that, every operation which is performed on target matrix, must also be performed to the identity matrix. In last, the target matrix will become identity matrix and the identity matrix will hold the inverse of the target matrix.

private static double determinant(double[,] matrix, int size)
    {
        double[] diviser = new double[size];// this will be used to make 0 all the elements of a row except (i,i)th value.
        double[] temp = new double[size]; // this will hold the modified ith row after divided by (i,i)th value.
        Boolean flag = false; // this will limit the operation to be performed only when loop n != loop i
        double determinant = 1;

        if (varifyRowsAndColumns(matrix, size)) // verifies that no rows or columns are similar or multiple of another row or column
            for (int i = 0; i < size; i++)
            {
                int count = 0;
                //this will hold the values to be multiplied by temp matrix
                double[] multiplier = new double[size - 1]; 
                diviser[i] = matrix[i, i];
                //if(i,i)th value is 0, determinant shall be 0
                if (diviser[i] == 0)
                {
                    determinant = 0;
                    break;
                }
                /*
                 * whole ith row will be divided by (i,i)th value and result will be stored in temp matrix.
                 * this will generate 1 at (i,i)th position in temp matrix i.e. ith row of matrix
                 */
                for (int j = 0; j < size; j++)
                {
                    temp[j] = matrix[i, j] / diviser[i];
                }
                //setting up multiplier to be used for multiplying the ith row of temp matrix
                for (int o = 0; o < size; o++)
                    if (o != i)
                        multiplier[count++] = matrix[o, i];

                count = 0;
                //for creating 0s at every other position than (i,i)th
                for (int n = 0; n < size; n++)
                {
                    for (int k = 0; k < size; k++)
                    {
                        if (n != i)
                        {
                            flag = true;
                            matrix[n, k] -= (temp[k] * multiplier[count]);
                        }
                    }
                    if (flag)
                        count++;
                    flag = false;
                }

            }
        else determinant = 0;
        //if determinant is not 0, (i,i)th element will be multiplied and the result will be determinant
        if (determinant != 0)
            for (int i = 0; i < size; i++)
            {
                determinant *= matrix[i, i];
            }
        return determinant;
    }


    private static Boolean varifyRowsAndColumns(double[,] matrix, int size)
    {
        List<double[]> rows = new List<double[]>();
        List<double[]> columns = new List<double[]>();


        for (int j = 0; j < size; j++)
        {
            double[] temp = new double[size];
            for (int k = 0; k < size; k++)
            {
                temp[j] = matrix[j, k];
            }
            rows.Add(temp);
        }

        for (int j = 0; j < size; j++)
        {
            double[] temp = new double[size];
            for (int k = 0; k < size; k++)
            {
                temp[j] = matrix[k, j];
            }
            columns.Add(temp);
        }

        if (!RowsAndColumnsComparison(rows, size))
            return false;
        if (!RowsAndColumnsComparison(columns, size))
            return false;
        return true;
    }
    private static Boolean RowsAndColumnsComparison(List<double[]> rows, int size)
    {
        int countEquals = 0;
        int countMod = 0;
        int countMod2 = 0;

        for (int i = 0; i < rows.Count; i++)
        {
            for (int j = 0; j < rows.Count; j++)
            {
                if (i != j)
                {
                    double min = returnMin(rows.ElementAt(i), rows.ElementAt(j));
                    double max = returnMax(rows.ElementAt(i), rows.ElementAt(j));

                    for (int l = 0; l < size; l++)
                    {
                        if (rows.ElementAt(i)[l] == rows.ElementAt(j)[l])
                            countEquals++;

                        for (int m = (int)min; m <= max; m++)
                        {
                            if (rows.ElementAt(i)[l] % m == 0 && rows.ElementAt(j)[l] % m == 0)
                                countMod++;

                            if (rows.ElementAt(j)[l] % m == 0 && rows.ElementAt(i)[l] % m == 0)
                                countMod2++;
                        }
                    }

                    if (countEquals == size)
                    {
                        return false;
                        // one row is equal to another row. determinant is zero
                    }

                    if (countMod == size)
                    {
                        return false;
                    }

                    if (countMod2 == size)
                    {
                        return false;
                    }
                }
            }
        }
        return true;
    }
    private static double returnMin(double[] row1, double[] row2)
    {
        double min1 = row1[0];
        double min2 = row2[0];

        for (int i = 1; i < row1.Length; i++)
            if (min1 > row1[i])
                min1 = row1[i];

        for (int i = 1; i < row2.Length; i++)
            if (min2 > row2[i])
                min2 = row2[i];

        if (min1 < min2)
            return min1;
        else return min2;
    }
    private static double returnMax(double[] col1, double[] col2)
    {
        double max1 = col1[0];
        double max2 = col2[0];

        for (int i = 1; i < col1.Length; i++)
            if (max1 < col1[i])
                max1 = col1[i];

        for (int i = 1; i < col2.Length; i++)
            if (max2 < col2[i])
                max2 = col2[i];

        if (max1 > max2)
            return max1;
        else return max2;
    }