Initial commit

diff --git a/matrix.c b/matrix.c
new file mode 100644 (file)
index 0000000..aeaa8f1
--- /dev/null
+++ b/matrix.c
@@ -0,0 +1,220 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include "matrix.h"
+
+void showMatrix(matrix *m){
+       for(int i = 0; i < m->rows; ++i){
+               for(int j = 0; j < m->cols; ++j){
+                       printf("%.1f\t", m->data[i][j]);
+               }
+               printf("\n");
+       }
+}
+
+matrix *newMatrix(int rows, int cols){
+       matrix *m = malloc(sizeof(matrix));
+       m->rows = rows;
+       m->cols = cols;
+       m->data = malloc(rows*sizeof(double*));
+       for(int i = 0; i < rows; ++i){
+               m->data[i] = malloc(cols*sizeof(double));
+       }
+       return m;
+}
+
+void freeMatrix(matrix **m){
+       for(int i = 0; i < (*m)->rows; ++i){
+               free((*m)->data[i]);
+               (*m)->data[i] = NULL;
+       }
+       free((*m)->data);
+       (*m)->data = NULL;
+       free(*m);
+       *m = NULL;
+}
+
+matrix *copyMatrix(matrix *m){
+       matrix *c = newMatrix(m->rows, m->cols);
+       for(int i = 0; i < m->rows; ++i){
+               for(int j = 0; j < m->cols; ++j){
+                       c->data[i][j] = m->data[i][j];
+               }
+       }
+       return c;
+}
+
+matrix *identityMatrix(int order){
+       matrix *im = newMatrix(order, order);
+       for(int i = 0; i < im->rows; ++i){
+               for(int j = 0; j < im->cols; ++j){
+                       im->data[i][j] = (i == j);
+               }
+       }
+       return im;
+}
+
+matrix *fillMatrix(matrix *m, double n){
+       matrix *r = newMatrix(m->rows, m->cols);
+       for(int i = 0; i < r->rows; ++i){
+               for(int j = 0; j < r->cols; ++j){
+                       r->data[i][j] = n;
+               }
+       }
+       return r;
+}
+
+matrix *addMatrix(matrix *m, double n){
+       for(int i = 0; i < m->rows; ++i){
+               for(int j = 0; j < m->cols; ++j){
+                       m->data[i][j] += n;
+               }
+       }
+}
+
+matrix *subtractMatrix(matrix *m, double n){
+       return addMatrix(m, -n);
+}
+
+matrix *scaleMatrix(matrix *m, double n){
+       matrix *r = copyMatrix(m);
+       for(int i = 0; i < r->rows; ++i){
+               for(int j = 0; j < r->cols; ++j){
+                       r->data[i][j] *= n;
+               }
+       }
+       return r;
+}
+
+static inline bool sameDimensions(matrix *m1, matrix *m2){
+       return (m1->rows == m2->rows) && (m1->cols == m2->cols);
+}
+
+matrix *addMatrices(matrix *m1, matrix *m2){
+       if(!sameDimensions(m1, m2)){
+               fprintf(stderr, "Wrong dimensions (%dx%d != %dx%d)\n", m1->rows, m1->cols, m2->rows, m2->cols);
+               return NULL;
+       }
+       matrix *r = newMatrix(m1->rows, m1->cols);
+       for(int i = 0; i < m1->rows; ++i){
+               for(int j = 0; j < m1->cols; ++j){
+                       r->data[i][j] = m1->data[i][j] + m2->data[i][j];
+               }
+       }
+       return r;
+}
+
+matrix *subtractMatrices(matrix *m1, matrix *m2){
+       m2 = scaleMatrix(m2, (double)-1);
+       return addMatrices(m1, m2);
+}
+
+matrix *multiplyMatrices(matrix *m1, matrix *m2){
+       if(m1->cols != m2->rows){
+               fprintf(stderr, "Wrong dimensions (%dx%d != %dx%d)\n", m1->rows, m1->cols, m2->rows, m2->cols);
+               return NULL;
+       }
+       matrix *r = newMatrix(m1->rows, m2->cols);
+       for(int i = 0; i < r->rows; ++i){
+               for(int j = 0; j < r->cols; ++j){
+                       double sum = 0;
+                       for(int n = 0; n < m1->cols; ++n){
+                               sum += m1->data[i][n] * m2->data[n][j];
+                       }
+                       r->data[i][j] = sum;
+               }
+       }
+       return r; 
+}
+
+static inline bool isSquare(matrix *m){
+       return m->rows == m->cols;
+}
+
+matrix *subMatrix(matrix *m, int row, int col){
+       matrix *r = newMatrix(m->rows-1, m->cols-1);
+       int ri = 0, rj = 0;
+       for(int i = 0; i < m->rows; ++i){
+               if(i == row){ continue; }
+               rj = 0;
+               for(int j = 0; j < m->cols; ++j){
+                       if(j == col){ continue; }
+                       r->data[ri][rj] = m->data[i][j];
+                       ++rj;
+               }
+               ++ri;
+       }
+       return r;
+}
+
+double determinant(matrix *m){
+       if(!isSquare(m)){
+               fprintf(stderr, "Matrix is not square (%dx%d)\n", m->rows, m->cols);
+               exit(EXIT_FAILURE);
+       }
+       double d = 0;
+       if(m->rows == 1){
+               return m->data[0][0];
+       }else{
+               for(int i = 0; i < m->rows; ++i){
+                       matrix *s = subMatrix(m, 0, i);
+                       double v = determinant(s)*m->data[0][i];
+                       d += !(i%2) ? v : -v;
+               }
+               return d;
+       }
+}
+
+matrix *cofactor(matrix *m){
+       if(!isSquare(m)){
+               fprintf(stderr, "Matrix is not square (%dx%d)\n", m->rows, m->cols);
+               return NULL;
+       }
+       matrix *r = newMatrix(m->rows, m->cols);
+       for(int i = 0; i < r->rows; ++i){
+               for(int j = 0; j < r->cols; ++j){
+                       double v = determinant(subMatrix(m, i, j));
+                       r->data[i][j] = !((i+j)%2) ? v : -v;
+               }
+       }
+       return r;
+}
+
+matrix *transpose(matrix *m){
+       matrix *r = newMatrix(m->cols, m->rows);
+       for(int i = 0; i < r->rows; ++i){
+               for(int j = 0; j < r->cols; ++j){
+                       r->data[i][j] = m->data[j][i];
+               }
+       }
+       return r;
+}
+
+matrix *adjugate(matrix *m){
+       return transpose(cofactor(m));
+}
+
+matrix *inverse(matrix *m){
+       double d = determinant(m);
+       if(d == 0){
+               fprintf(stderr, "Determinant is 0, the matrix is not invertible\n");
+               return NULL;
+       }
+       return scaleMatrix(adjugate(m), 1/d);
+}
+
+matrix *raiseMatrix(matrix *m, int n){
+       if(!isSquare(m)){
+               fprintf(stderr, "Matrix is not square (%dx%d)\n", m->rows, m->cols);
+               return NULL;
+       }
+       if(n < 0){
+               m = inverse(m);
+               n = -n;
+       }
+       matrix *r = identityMatrix(m->rows);
+       for(int i = 0; i < n; ++i){
+               r = multiplyMatrices(r, m);
+       }
+       return r;
+}

diff --git a/matrix.h b/matrix.h
new file mode 100644 (file)
index 0000000..b01c889
--- /dev/null
+++ b/matrix.h
@@ -0,0 +1,48 @@
+#pragma once
+#ifndef MATRIX_H
+#define MATRIX_H
+
+
+typedef struct matrix{
+       int rows;
+       int cols;
+       double **data;
+} matrix;
+
+matrix *newMatrix(int rows, int cols);
+
+void freeMatrix(matrix **m);
+
+matrix *copyMatrix(matrix *m);
+
+matrix *identityMatrix(int order);
+
+matrix *fillMatrix(matrix *m, double n);
+
+matrix *addMatrix(matrix *m, double n);
+
+matrix *subtractMatrix(matrix *m, double n);
+
+matrix *scaleMatrix(matrix *m, double n);
+
+matrix *addMatrices(matrix *m1, matrix *m2);
+
+matrix *subtractMatrices(matrix *m1, matrix *m2);
+
+matrix *multiplyMatrices(matrix *m1, matrix *m2);
+
+matrix *subMatrix(matrix *m, int row, int col);
+
+double determinant(matrix *m);
+
+matrix *cofactor(matrix *m);
+
+matrix *transpose(matrix *m);
+
+matrix *adjugate(matrix *m);
+
+matrix *inverse(matrix *m);
+
+matrix *raiseMatrix(matrix *m, int n);
+
+#endif