#include <stdio.h>
#include <math.h>

// Prototypes
void mxinverse(double* matrix, int n);
void mxmult(double* matrixA, double* matrixB, double* matrixC,
            int mA, int nA, int mB, int nB);
void mxadd(double* matrixA, double* matrixB, double* matrixC,
            int mA, int nA, int mB, int nB, double alpha, double beta);


// ToDo: Improve this function returning a status value
void mxinverse(double* matrix, int n) {
    int i, j, k;
    double sum;

    for (i = 0; i < n; i++) {
        for (j = i; j < n; j++) {
            sum = matrix[i * n + j];
            for (k = i - 1; k >= 0; k--) {
                sum -= matrix[i * n + k] * matrix[j * n + k];
            }
            if (i == j) {
                // What happen when sum < 0.0?
                matrix[i * n + i] = sqrt(sum);
            } else {
                // What happen when matrix[i * n + i] == 0.0?
                // ToDo: return a status value
                matrix[j * n + i] = sum / matrix[i * n + i];
            }
        }
    }

    for (i = 0; i < n; i++) {
        for (j = i; j < n; j++) {
            sum = (i == j) ? 1.0 : 0.0;
            for (k = j - 1; k >= i; k--) {
                sum -= matrix[j * n + k] * matrix[j * n + k];
            }
            // What happen when matrix[i * n + i] == 0.0?
            // ToDo: return a status value
            matrix[j * n + i] = matrix[j * n + i] / matrix[i * n + i];
            matrix[i * n + j] = sum / matrix[i * n + i];
        }
    }

    for (i = 0; i < n; i++) {
        for (j = i + 1; j < n; j++) {
            for (k = i + 1; k < j; k++) {
                matrix[j * n + i] -= matrix[j * n + k] * matrix[k * n + i];
            }
                // What happen when matrix[i * n + i] == 0.0?
                // ToDo: return a status value
                matrix[j * n + i] /= matrix[j * n + j];
        }
    }
}


// C[i,j] = alpha * ( A[i,j] * B[i,j] ) + beta * C[i,j] // C = A  * B
// C[i,j] = alpha * ( A[j,i] * B[i,j] ) + beta * C[i,j] // C = A' * B
// C[i,j] = alpha * ( A[j,i] * B[j,i] ) + beta * C[i,j] // C = A' * B'
// C[i,j] = alpha * ( A[i,j] * B[j,i] ) + beta * C[i,j] // C = A  * B'
// 
void mxmult(double* matrixA, double* matrixB, double* matrixC,
            int mA, int nA, int mB, int nB, double alpha, double beta) {               
  // ToDo: Implement the algorithm
}


// C[i,j] = alpha * A[i,j] + beta * B[i,j] + gamma
//
// Changing alpha and beta, it is possible to implement mxsub
void mxadd(double* matrixA, double* matrixB,
           int mA, int nA, int mB, int nB,
           double alpha, double beta, double gamma) {
    
}