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

#include "mxlib.h"

/** Please note that the indexing of rows and columns always starts at 0 */
/* For example, the data in the upper left corner of a matrix(M) is addressed as M->d[0][0] */

/** Every function in the library requires you to pass matrices by address to them, and will return pointers to matrices (if applicable) */

/** The Matrix_adv_ matrix operation functions can be used to delete matrices passed to them before returning */
/* An example where this can be used: */
/* Matrix *C = Matrix_adv_add(1, Matrix_transpose(A), 0, B); */
/* Here, without using the _adv_ version of the add function, the pointer to the transposed A matrix would be thrown away, thus causing memory leak */
/* These functions are separate from their simple counterparts to reduce unnecessary clutter of 0's and 1's in every matrix function call */
/* Another example can be found in the included main.c file */

/** There is an included main.c file with some examples */

/** All non-trivial functions are commented in this (mxlib.c) source file */

/** All error messages are printed to stderr */
/* Error messages start with "mxlib:" */
/* I implemented error checking wherever I could think to do so, I hope I covered most problems that can occur */

/** Please contact developer (zsmb13@sch.bme.hu) with any feedback */


/** ------------------------------------------------- */
/** --------- STATIC FUNCTION DECLARATIONS ---------- */
/** ------------------------------------------------- */

static Matrix *Matrix_allocate(int k, int n);
static int Matrix_intdigits(Matrix *M);
static void double_switch(double *a, double *b);
static int switch_possible(Matrix *M, int k, int n, int i, int j);
static int non_zero_constant(Matrix *M, int k, int n, int i);
static Matrix *Matrix_Gaussian_elimination_phase_1(Matrix *M, int k, int n, int *solvable, int lin);
static Matrix *Matrix_Gaussian_elimination_phase_2(Matrix *M, int k, int n);
static Matrix *Matrix_Gaussian_elimination(Matrix *M, int k, int n, int lin, int *solutions);
static void Matrix_column_copy(Matrix *A, int j1, Matrix *B, int j2);
static Matrix *Matrix_half(Matrix *M);

/** ------------------------------------------------- */
/** ---------------- MATRIX CREATION ---------------- */
/** ------------------------------------------------- */

/* Function used by the Matrix_create_ functions for memory allocation */
static Matrix *Matrix_allocate(int k, int n) {
    if(k <= 0 || n <= 0) {
        fprintf(stderr, "mxlib: The minimum size of a matrix is 1x1!\n");
        return NULL;
    }

    Matrix *M = malloc(sizeof(Matrix));
    if(M == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for a matrix has failed!\n");
        return NULL;
    }
    int i;

    M->k = k;
    M->n = n;
    M->d = malloc(k*sizeof(double*));
    for(i = 0; i < k; i++)
        M->d[i] = malloc(n*sizeof(double));

    return M;
}

/* Frees all memory allocated to a matrix */
Matrix *Matrix_delete(Matrix *M) {
    if(M == NULL) {
        fprintf(stderr, "mxlib: NULL matrix passed to the Matrix_delete function!\n");
        return NULL;
    }
    int i;
    for(i = 0; i < M->k; i++)
        free(M->d[i]);
    free(M->d);
    free(M);
    return NULL; /* Return value can be used for setting pointer to NULL */
}

/* Creates a matrix of size kxn with all data points set to 0 */
Matrix *Matrix_create_nullmatrix(int k, int n) {
    if(k <= 0 || n <= 0) {
        fprintf(stderr, "mxlib: The minimum size of a matrix is 1x1!\n");
        return NULL;
    }

    Matrix *M = Matrix_allocate(k, n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Null matrix creation failed!\n");
        return NULL;
    }

    int i, j;

    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++)
            M->d[i][j] = 0;

    return M;
}

/* Creates a unit matrix of size nxn */
Matrix *Matrix_create_unitmatrix(int n) {
    if(n <= 0) {
        fprintf(stderr, "mxlib: The minimum size of a matrix is 1x1!\n");
        return NULL;
    }

    Matrix *M = Matrix_create_nullmatrix(n, n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Unit matrix creation failed!\n");
        return NULL;
    }

    int i;

    for(i = 0; i < n; i++)
        M->d[i][i] = 1;

    return M;
}

/* All data entered after the first two arguments must be of type double, cast integers to double when calling this function, or add .0 when manually creating a matrix! */
Matrix *Matrix_create(int k, int n, ...) {
    if(k <= 0 || n <= 0) {
        fprintf(stderr, "mxlib: The minimum size of a matrix is 1x1!\n");
        return NULL;
    }

    double data;
    va_list dp;
    int i = 0, j = 0;
    Matrix *M = Matrix_allocate(k, n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_create has failed!\n");
        return NULL;
    }

    va_start(dp, n);

    while(i < k) {
        data = va_arg(dp, double);
        M->d[i][j] = data;
        j++;
        if(j == n) {
            i++;
            j = 0;
        }
    }

    va_end(dp);
    return M;
}

/** ------------------------------------------------- */
/** -------------- DISPLAYING MATRICES -------------- */
/** ------------------------------------------------- */

/* This function is used in the Matrix_display function to determine the number of integer digits it has to display */
static int Matrix_intdigits(Matrix *M) {
    int i, j;
    int roundedvalue;
    int digits = 0;
    double maxvalue = fabs(M->d[0][0]);

    /* Looks up the data point with the largest absolute value in the matrix */
    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++) {
            if(fabs(M->d[i][j]) > maxvalue)
                maxvalue = fabs(M->d[i][j]);
        }
    roundedvalue = (int) ceil(maxvalue);

    /* Calculates number of digits */
    while(roundedvalue) {
        roundedvalue /= 10;
        digits++;
    }

    return digits;
}

/* Precision variable determines number of fractional digits displayed */
/* Number of integer digits displayed is calculated automatically */
void Matrix_display(Matrix *M, int precision) {
    if(M == NULL) {
        fprintf(stderr, "mxlib: NULL matrix passed to the Matrix_display function!\n");
        return;
    }
    if(precision < 0) {
        fprintf(stderr, "mxlib: The precision variable has to be a non-negative integer!\n");
        return;
    }

    int i, j;
    int digits = 1 + Matrix_intdigits(M) + ((precision>0)?1:0) + precision; /* +/- sign, integer digits, decimal point (1 char), fractional digits */
    for(i = 0; i < M->k; i++) {
        for(j = 0; j < M->n; j++)
            printf("% *.*f ", digits, precision, (fabs(M->d[i][j]) < pow(10, -13))?0.0:M->d[i][j]); /* Displays extremely small absolute value numbers as 0 */
        printf("\n");
    }
}

/* Displays matrix with a separator line before the column that's index is passed in as the last parameter */
/* Can be used for displaying augmented matrices */
void Matrix_display_with_separator(Matrix *M, int precision, int separator_placement) {
    if(M == NULL) {
        fprintf(stderr, "mxlib: NULL matrix passed to the Matrix_display function!\n");
        return;
    }
    if(precision < 0) {
        fprintf(stderr, "mxlib: The precision variable has to be a non-negative integer!\n");
        return;
    }

    int i, j;
    int digits = 1 + Matrix_intdigits(M) + ((precision>0)?1:0) + precision; /* +/- sign, integer digits, decimal point (1 char), fractional digits */
    for(i = 0; i < M->k; i++) {
        for(j = 0; j < M->n; j++) {
            if(j == separator_placement) {
                putchar('|');
                putchar(' ');
            }
            printf("% *.*f ", digits, precision, (fabs(M->d[i][j]) < pow(10, -13))?0.0:M->d[i][j]); /* Displays extremely small absolute value numbers as 0 */
        }
        printf("\n");
    }
}

/** ------------------------------------------------- */
/** ---------------- FILE OPERATIONS ---------------- */
/** ------------------------------------------------- */

/* Writes the values of a matrix into a text file */
void Matrix_file_write(char *filename, Matrix *M) {
    if(M == NULL) {
        fprintf(stderr, "mxlib: NULL matrix passed to the Matrix_file_write function!\n");
        return;
    }

    FILE *fp = fopen(filename, "wt");
    if(fp == NULL) {
        fprintf(stderr, "mxlib: Failed to open file for writing in function Matrix_file_write !\n");
    }

    int i, j;
    for(i = 0; i < M->k; i++) {
        for(j = 0; j < M->n; j++)
            fprintf(fp, "%f ", M->d[i][j]);
        fprintf(fp, "\n");
    }

    fclose(fp);
}

/* Builds a matrix from values in a text file */
/* Can not do error checking for the k and n values! */
Matrix *Matrix_file_read(char *filename, int k, int n) {
    FILE *fp = fopen(filename, "rt");
    if(fp == NULL) {
        fprintf(stderr, "mxlib: Failed to open file for reading in function Matrix_file_read !\n");
    }

    Matrix *M = Matrix_allocate(k, n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_file_read has failed!\n");
        return NULL;
    }

    int i, j;
    for(i = 0; i < M->k; i++) {
        for(j = 0; j < M->n; j++)
            fscanf(fp, "%lf ", &M->d[i][j]);
        //fscanf(fp, "\n");
    }

    fclose(fp);

    return M;
}

/** ------------------------------------------------- */
/** ------------- MATRIX ROW OPERATIONS ------------- */
/** ------------------------------------------------- */
/* All of these functions overwrite the original matrix */

/* Multiply all values in a row of a matrix by the constant c */
Matrix *Matrix_row_multiply(Matrix *M, int i, double c) {
    if(i < 0 || M->k < i) {
        fprintf(stderr, "mxlib: No such row to multiply!\n");
        return M;
    }

    int j;
    for(j = 0; j < M->n; j++)
        M->d[i][j] = M->d[i][j] * c;
    return M;
}

/* Adds a row of a matrix (i2) to another (i1) */
Matrix *Matrix_row_add(Matrix *M, int i1, int i2) {
    if(i1 < 0 || M->k < i1 || i2 < 0 || M->k < i2) {
        fprintf(stderr, "mxlib: No such rows to add!\n");
        return M;
    }

    int j;
    for(j = 0; j < M->n; j++)
        M->d[i1][j] += M->d[i2][j];

    return M;
}

/* Multiples a row of a matrix (i2) by the constant c, then adds it to another (i1) */
Matrix *Matrix_row_multiply_and_add(Matrix *M, int i1, int i2, double c) {
    if(i1 < 0 || M->k < i1 || i2 < 0 || M->k < i2) {
        fprintf(stderr, "mxlib: No such rows to add and multiply!\n");
        return M;
    }

    int j;
    for(j = 0; j < M->n; j++)
        M->d[i1][j] += M->d[i2][j] * c;

    return M;
}

/* Helper function for the Matrix_row_switch function */
/* Swaps the values of two doubles passed by address */
static void double_switch(double *a, double *b) {
    double temp = *a;
    *a = *b;
    *b = temp;
}

/* Switches up two rows of a matrix */
Matrix *Matrix_row_switch(Matrix *M, int i1, int i2) {
    if(i1 < 0 || M->k < i1 || i2 < 0 || M->k < i2) {
        fprintf(stderr, "mxlib: No such rows to switch!\n");
        return M;
    }

    int j;
    for(j = 0; j < M->n; j++)
        double_switch(&M->d[i1][j], &M->d[i2][j]);

    return M;
}

/* Deletes a row from a matrix */
Matrix *Matrix_row_delete(Matrix *M, int i) {
    if(i < 0 || M->k < i) {
        fprintf(stderr, "mxlib: No such row to delete!\n");
        return M;
    }

    for(; i < M->k-1; i++) {
        M = Matrix_row_switch(M, i, i+1);
    }

    free(M->d[i]);

    M->k = M->k - 1;
    return M;
}

/** ------------------------------------------------- */
/** -------------- GAUSSIAN ELIMINATION ------------- */
/** ------------------------------------------------- */

/* Short functions for checking things during the elimniation*/
/* Looks for a row that has a non-zero value below the current position */
static int switch_possible(Matrix *M, int k, int n, int i, int j) {
    int t;
    for(t = i+1; t <= k; t++)
        if(M->d[t][j] != 0)
            return t;
    return 0;
}

/* Checks if there are non-zero constants in all-zero rows */
/* Used only if the elimination is ran to solve a system of linear equations */
static int non_zero_constant(Matrix *M, int k, int n, int i) {
    int t;
    for(t = i+1; t <= k; t++)
        if(M->d[t][n+1] != 0)
            return 1;
    return 0;
}

/* Algorithm used is from http://cs.bme.hu/bsz1/jegyzet/bsz1_jegyzet.pdf , pages 37-39 */
/* First phase of the elimination, reaching the low echelon form */
static Matrix *Matrix_Gaussian_elimination_phase_1(Matrix *M, int k, int n, int *solvable, int lin) {
    int t;
    int i = 0, j = 0;

    step1: /* Normal procedure */
        if(M->d[i][j] == 0) goto step2;

        M = Matrix_row_multiply(M, i, 1/M->d[i][j]);

        if(i == k) goto step3;

        for(t = i+1; t <= k; t++)
            M = Matrix_row_multiply_and_add(M, t, i, (-1)*M->d[t][j]);

        if(j == n) goto step3;

        i++, j++;

        goto step1;

    step2: /* Exception handling */
        if(i < k && (t = switch_possible(M, k, n, i, j))) {
            M = Matrix_row_switch(M, i, t); goto step1;}

        if(j == n) {
            i--; goto step3;}

        j++;

        goto step1;

    step3: /* Done, final steps */
        if(i == k) {*solvable = 1; return M;}

        if(lin && non_zero_constant(M, k, n, i)) {*solvable = 0; return M;}

        for(t = i+1; t <= k; t++)
            M = Matrix_row_delete(M, i+1);

         *solvable = 2; return M;
}

/* Conducts the second phase of the elimination, computes the matrix to a reduced row echelon form */
static Matrix *Matrix_Gaussian_elimination_phase_2(Matrix *M, int k, int n) {
    int t;
    int i, j;

    for(i = k; i > 0; i--) {
        j = 0;
        while(M->d[i][j] != 1)
            j++;
        for(t = i-1; t >= 0; t--)
            M = Matrix_row_multiply_and_add(M, t, i, -1*M->d[t][j]);
    }

    return M;
}

/* The values of k and n determine the size of the matrix the elimination should be computed on */
/* This function doesn't get called by the user, only by other functions that know what k and n parameters to use */
/* Should be called with (0, NULL) as last parameters in all functions other than the linear equation solver */
static Matrix *Matrix_Gaussian_elimination(Matrix *M, int k, int n, int lin, int *solutions) {
    int solvable;
    k = k - 1;
    n = n - 1;

    M = Matrix_Gaussian_elimination_phase_1(M, k, n, &solvable, lin);

    switch(solvable) {
        case 0:
            //printf("No solutions");
        break;
        case 1:
        case 2:
            M = Matrix_Gaussian_elimination_phase_2(M, M->k-1, M->n-1);
        break;
    }

    if(solutions != NULL)
        *solutions = solvable;

    return M;
}

/* Solves a system of linear equations using Gaussian elimination */
/* Overwrites original with the reduced low echelon form matrix */
/* Return values: */
/* 0 - no solutions */
/* 1 - one solution */
/* 2 - infinite solutions */
int Matrix_linear_equations_solve(Matrix *M) {
    int solutions;
    Matrix_Gaussian_elimination(M, M->k, M->n-1, 1, &solutions);
    return solutions;
}


/* Algorithm used is from http://cs.bme.hu/bsz1/jegyzet/bsz1_jegyzet.pdf , page 52 */
/* Included in this section because it uses a simpler version of Gaussian elimination */
double Matrix_determinant(Matrix *M) {
    if(M->k != M->n) {
        fprintf(stderr, "mxlib: Determinant is only defined for square matrices!\n");
        return 0;
    }

    /* If matrix is upper or lower triangular, the determinant is the product of the values in the main diagonal of the matrix */
    if(Matrix_is_lower_triangular(M) || Matrix_is_upper_triangular(M)) {
        int i, det = 1;
        for(i = 0; i < M->k; i++)
            det *= M->d[i][i];
        return det;
    }

    int t;
    double det = 1;
    int i = 0;
    int n = M->n - 1;
    Matrix *C = Matrix_copy(M);

    step1:
        if(C->d[i][i] == 0) goto step2;

        det *= C->d[i][i];

        C = Matrix_row_multiply(C, i, 1/C->d[i][i]);

        if(i == n) {Matrix_delete(C); return det;}

        for(t = i+1; t <= n; t++)
            C = Matrix_row_multiply_and_add(C, t, i, (-1)*C->d[t][i]);

        i++;

        goto step1;

    step2:
        if(i < n && (t = switch_possible(C, n, n, i, i))) {
           C = Matrix_row_switch(C, i, t);
           det *= (-1);
           goto step1;
        }
        Matrix_delete(C);
        return 0.0;
}

/** ------------------------------------------------- */
/** -------------- MISC MATRIX FUNCTIONS ------------ */
/** ------------------------------------------------- */

/* Determines whether a matrix is an upper triangular matrix */
int Matrix_is_upper_triangular(Matrix *M) {
    if(M->k != M->n)
        return 0;
    int i, j;
    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++)
            if(i > j && M->d[i][j] != 0)
                return 0;
    return 1;
}

/* Determines whether a matrix is a lower triangular matrix */
int Matrix_is_lower_triangular(Matrix *M) {
    if(M->k != M->n)
        return 0;
    int i, j;
    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++)
            if(i < j && M->d[i][j] != 0)
                return 0;
    return 1;
}

/* Returns the rank of the matrix, determined by counting the number of rows left after a Gaussian elimination */
/* Doesn't modify the original matrix in the process */
int Matrix_rank(Matrix *M) {
    int r;
    Matrix *C = Matrix_copy(M);
    C = Matrix_Gaussian_elimination(C, C->k, C->n, 0, NULL);
    r = C->k;
    Matrix_delete(C);
    return r;
}

/* Used in the Matrix_augment functions */
static void Matrix_column_copy(Matrix *A, int j1, Matrix *B, int j2) {
    if(A->k != B->k) {
        fprintf(stderr, "mxlib: Column copy failed, number of rows mismatch!\n");
        return;
    }

    int i;
    for(i = 0; i < A->k; i++)
        A->d[i][j1] = B->d[i][j2];
}

/* Used in the Matrix_inverse function (gets the right half of the augmented matrix) */
static Matrix *Matrix_half(Matrix *M) {
    Matrix *H = Matrix_allocate(M->k, M->n/2);
    int i;
    for(i = 0; i < H->n; i++)
        Matrix_column_copy(H, i, M, i+H->n);
    Matrix_delete(M);
    return H;
}

/** ------------------------------------------------- */
/** --------------- MATRIX OPERATIONS --------------- */
/** ------------------------------------------------- */

/* Returns a freshly allocated copy of a matrix */
Matrix *Matrix_copy(Matrix *M) {
    if(M == NULL) return NULL;

    int i, j;
    Matrix *C = Matrix_allocate(M->k, M->n);
    if(C == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_copy has failed!\n");
        return NULL;
    }

    C->k = M->k;
    C->n = M->n;
    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++)
            C->d[i][j] = M->d[i][j];

    return C;
}

Matrix *Matrix_add(Matrix *A, Matrix *B) {
    if(A->k != B->k || A->n != B->n) {
        fprintf(stderr, "mxlib: Matrix addition is only defined for two matrices of the same size!\n");
        return NULL;
    }

    int i, j;
    Matrix *M = Matrix_allocate(A->k, A->n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_add has failed!\n");
        return NULL;
    }
    M->k = A->k;
    M->n = A->n;
    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++)
            M->d[i][j] = A->d[i][j] + B->d[i][j];

    return M;
}

Matrix *Matrix_subtract(Matrix *A, Matrix *B) {
    if(A->k != B->k || A->n != B->n) {
        fprintf(stderr, "mxlib: Matrix subtraction is only defined for two matrices of the same size!\n");
        return NULL;
    }

    int i, j;
    Matrix *M = Matrix_allocate(A->k, A->n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_subtract has failed!\n");
        return NULL;
    }
    M->k = A->k;
    M->n = A->n;
    for(i = 0; i < M->k; i++)
        for(j = 0; j < M->n; j++)
            M->d[i][j] = A->d[i][j] - B->d[i][j];

    return M;
}

Matrix *Matrix_multiply_by_scalar(Matrix *M, double c) {
    if(M == NULL) {
        fprintf(stderr, "mxlib: NULL matrix passed to the Matrix_multiply_by_scalar function!\n");
        return NULL;
    }
    int i, j;
    Matrix *cM = Matrix_allocate(M->k, M->n);
    if(cM == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_multiply_by_scalar has failed!\n");
        return NULL;
    }
    cM->k = M->k;
    cM->n = M->n;
    for(i = 0; i < cM->k; i++)
        for(j = 0; j < cM->n; j++)
            cM->d[i][j] = c * M->d[i][j];

    return cM;
}

Matrix *Matrix_transpose(Matrix *M) {
    if(M == NULL) {
        fprintf(stderr, "mxlib: NULL matrix passed to the Matrix_transpose function!\n");
        return NULL;
    }

    int i, j;
    Matrix *MT = Matrix_allocate(M->n, M->k);
    if(MT == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_transpose has failed!\n");
        return NULL;
    }
    MT->k = M->n;
    MT->n = M->k;
    for(i = 0; i < MT->k; i++)
        for(j = 0; j < MT->n; j++)
            MT->d[i][j] = M->d[j][i];

    return MT;
}

Matrix *Matrix_multiply(Matrix *A, Matrix *B) {
    if(A->n != B->k) {
        fprintf(stderr, "mxlib: Dimension error when calling the Matrix_multiply function!\n");
        fprintf(stderr, "mxlib: (Number of columns in A does not match number of rows in B)\n");
        return NULL;
    }

    int h, i, j;
    Matrix *AB = Matrix_create_nullmatrix(A->k, B->n);
    if(AB == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_multiply has failed!\n");
        return NULL;
    }

    for(i = 0; i < AB->k; i++)
        for(h = 0; h < A->n; h++)
            for(j = 0; j < AB->n; j++)
                AB->d[i][j] += A->d[i][h] * B->d[h][j];
    return AB;
}

Matrix *Matrix_inverse(Matrix *M) {
    if(M->k != M->n) {
        fprintf(stderr, "mxlib: Only square matrices have inverses!\n");
        return M;
    }
    Matrix *I = Matrix_adv_augment(0, M, 1, Matrix_create_unitmatrix(M->n));
    I = Matrix_Gaussian_elimination(I, I->k, I->n/2, 0, NULL);
    I = Matrix_half(I);
    return I;
}

/* Augments two matrices to each other*/
Matrix *Matrix_augment(Matrix *A, Matrix *B) {
    if(A->k != B->k) {
        fprintf(stderr, "mxlib: Augmentation failed, number of rows mismatch!\n");
        return NULL;
    }
    Matrix *M = Matrix_allocate(A->k, A->n + B->n);
    if(M == NULL) {
        fprintf(stderr, "mxlib: Memory allocation for function Matrix_augment has failed!\n");
        return NULL;
    }

    int i;
    for(i = 0; i < A->k; i++)
        Matrix_column_copy(M, i, A, i);
    for(; i < A->n + B->n; i++)
        Matrix_column_copy(M, i, B, i-A->k);

    return M;
}

/** ------------------------------------------------- */
/** ----------- ADVANCED MATRIX FUNCTIONS ----------- */
/** ------------------------------------------------- */

void Matrix_adv_display(int delM, Matrix *M, int precision) {
    Matrix_display(M, precision);
    if(delM) Matrix_delete(M);
}

void Matrix_adv_display_with_separator(int delM, Matrix *M, int precision, int separator_placement) {
    Matrix_display_with_separator(M, precision, separator_placement);
    if(delM) Matrix_delete(M);
}

Matrix *Matrix_adv_add(short delA, Matrix *A, short delB, Matrix *B) {
    Matrix *M = Matrix_add(A, B);
    if(delA) Matrix_delete(A);
    if(delB) Matrix_delete(B);
    return M;
}

Matrix *Matrix_adv_subtract(short delA, Matrix *A, short delB, Matrix *B) {
    Matrix *M = Matrix_subtract(A, B);
    if(delA) Matrix_delete(A);
    if(delB) Matrix_delete(B);
    return M;
}

Matrix *Matrix_adv_multiply_by_scalar(short delM, Matrix *M, double c) {
    Matrix *cM = Matrix_multiply_by_scalar(M, c);
    if(delM) Matrix_delete(M);
    return cM;
}

Matrix *Matrix_adv_transpose(short delM, Matrix *M) {
    Matrix *MT = Matrix_transpose(M);
    if(delM) Matrix_delete(M);
    return MT;
}

Matrix *Matrix_adv_multiply(short delA, Matrix *A, short delB, Matrix *B) {
    Matrix *AB = Matrix_multiply(A, B);
    if(delA) Matrix_delete(A);
    if(delB) Matrix_delete(B);
    return AB;
}

Matrix *Matrix_adv_inverse(int delM, Matrix *M) {
    Matrix *I = Matrix_inverse(M);
    if(delM) Matrix_delete(M);
    return I;
}

Matrix *Matrix_adv_augment(int delA, Matrix *A, int delB, Matrix *B) {
    Matrix *M = Matrix_augment(A, B);
    if(delA) Matrix_delete(A);
    if(delB) Matrix_delete(B);
    return M;
}