Thread: Matrice Recursive Determinant Calculation

maybe like this (not tested for bigger matrix yet, sorry)

Code:

#include <stdio.h>

#define MY_FREE(x)  { free(x); x = NULL; }

typedef struct {
    unsigned int width;
    unsigned int height;
    float        *data;
} matrix;

/*
 * new_matrix : creates new matrix
 * w          : width (column)
 * h          : height (row)
 * returns    : newly allocated matrix (on success), NULL (otherwise)
 */
matrix *new_matrix (unsigned int w, unsigned int h)
{
    matrix *result = (matrix *)malloc(sizeof(matrix));

    if (result != NULL) {
	result->width  = w;
	result->height = h;
	result->data   = (float *)calloc(w * h, sizeof(float));
	if (result->data == NULL) {
	    fprintf(stderr, "new_matrix() : out of memory.\n");
	    MY_FREE(result);
	}
    }

    return(result);
}

/*
 * destroy_matrix : it's obvious
 */
void destroy_matrix (matrix *mat)
{
    if (mat != NULL) {
	MY_FREE(mat->data);
	MY_FREE(mat);
    }
}

/* 
 * determinant : get the determinant of a matrix
 * mat     : matrix input
 * retval  : 1 (success), 0 (failed)
 * returns : determinant of matrix mat
 */
float determinant (const matrix *mat, int *retval)
{
    float result;

    if (mat->width != mat->height) {
	fprintf(stderr, "not a square matrix!\n");
	*retval = 0;	/* return false */
    } else {
	if (mat->width == 2) {
	    /* det = ad - bc */
	    result  = (mat->data[0] * mat->data[3]) - (mat->data[1] * mat->data[2]);
	    *retval = 1;
	} else {
	    matrix *tmp_mat = new_matrix(mat->width - 1, mat->height - 1);

	    if (tmp_mat == NULL) {
		*retval = 0;
	    } else {
		const float sign[2] = { 1.0, -1.0 };
		int         col, rc, x, y, xx, yy;
		float       tmp_det;

		result = 0.0;
		for (col = 0; col < mat->width; ++col) {
		    for (y = 1, yy = 0; y < mat->height; ++y, ++yy) {
			for (x = 0, xx = 0; x < mat->width; ++x) {
			    if (x != col) {
				tmp_mat->data[yy * tmp_mat->width + xx] = mat->data[y * mat->width + x];
				++xx;
			    }
			}
		    }
		    tmp_det = determinant(tmp_mat, &rc);
		    if (rc) {	/* success */
			result = result + (mat->data[col] * sign[col &#37; 2] * tmp_det);
		    } else {	/* failed */
			*retval = 0;
			destroy_matrix(tmp_mat);
			break;
		    }
		}
		*retval = 1;
	    }
	}
    }

    return result;
}

/* test program */
int main (void)
{
    matrix *mat1 = new_matrix(3, 3);
    int    x, y, rc;
    float  det;

    for (x = 0; x < mat1->width * mat1->height; ++x) {
	mat1->data[x] = (float)1;
    }
    mat1->data[4] = 5.0;
    mat1->data[6] = 13.0;
    det = determinant(mat1, &rc);
    if (rc) {
	printf("determinant : %f\n", det);
    }
    for (y = 0; y < mat1->width; ++y) {
	for (x = 0; x < mat1->height; ++x) {
	    printf("%5.2f", mat1->data[y * mat1->width + x]);
	}
	printf("\n");
    }
    destroy_matrix(mat1);

    return(0);
}