determinant algorithm of a 4x4 matrix

Question

I pick the first row and multiply each element by its cofactor, but in some cases the method is returning nan. For example,

1 0 0 1
0 2 0 0
0 0 3 0
0 0 0 4

in this case the method returns nan.

Does anyone know what I did wrong?

getDet3 returns determinant of a 3x3 matrix and it works fine.

-(double) getDet4:(double[4][4])mat {
    double det = 0;
    double small[3][3];

    int i, j, k;
    int i_ = 1, j_;

    for ( i=0; i<4; i++ ){

        if (mat[0][i] == 0) continue;

        // get the small matrix here
        for ( j=0; j<3; j++ ){
            j_ = 0;
            for ( k=0; k<3; k++ ){
                if ( i == j_ ) j_++;
                small[j][k] = mat[i_][j_];
                j_++;
            }
            i_++;
        }

        det += mat[0][i] * [self getDet3:small] * pow(-1, i+j);
    }

    return det;
}

Answer 1

Well, there are a few mistakes in your code.

1) The initialization of i_ = 1 should be done just before the j loop, otherwise it will keep the old value.

2) The computation of pow(-1, i+j) should only depend on i, since j has the same value every time in that expression (namely, 3).

So, assuming that getDet3 is correct, the mistake is introduced by i_ going out of bounds. As a whole, the code should look like:

-(double) getDet4:(double[4][4])mat {
    double det = 0;
    double small[3][3];

    int i, j, k;
    int i_, j_;

    for ( i=0; i<4; i++ ){

        if (mat[0][i] == 0) continue;

        // get the small matrix here
        i_ = 1;
        for ( j=0; j<3; j++ ){
            j_ = 0;
            for ( k=0; k<3; k++ ){
                if ( i == j_ ) j_++;
                small[j][k] = mat[i_][j_];
                j_++;
            }
            i_++;
        }

        det += mat[0][i] * [self getDet3:small] * pow(-1, i);
    }

    return det;
}

Answer 2

Personally, I find your variable names confusing. If I understand your idea correctly, you expect i_ to have the value j + 1 and j_ to be k < i ? k : k + 1. IMHO, it would have been less confusing to have named them j_p andk_`, or even to just use the equivalent expression.

In any event, you don't reinitialize i_ inside the outer for loop. So it actually just keeps on incrementing, resulting in array indices outside of the array bounds.