#ifndef MAT_H
#define MAT_H

#include <gmp.h>
#include <malloc.h>

#define LOOP(i,n) for(int i = 0; i < (n); i++)

// library for matrix computations in variable rings (based on GMP types)

/*
  needed features:
  x multiply matrices
  - inverse
  - pseudoinverse
  x set
  - eigenvalues
 */

#define NUMBER mpq_t
#define INIT mpq_init
#define CLEAR mpq_clear
#define SET mpq_set
#define SET_ZERO(x) mpq_set_ui(x,0,1)
#define SET_ONE(x) mpq_set_ui(x,1,1)
#define ADD mpq_add
#define SUB mpq_sub
#define MULTIPLY mpq_mul
#define DIV mpq_div
#define PRINT(x) gmp_printf("%Qd", x)

#define M(m,i,j) ((m)->x[(i)+(m)->n*(j)])

struct _mat{
	int n;
	NUMBER *x;
} ;

typedef struct _mat mat[1];

typedef struct _mat_workspace {
	mat tmp_mat;
	NUMBER tmp_num;
	NUMBER tmp_num2;
} mat_workspace;

mat_workspace *mat_workspace_init(int n);
void mat_workspace_clear(mat_workspace *ws);
void mat_init(mat m, int n);
void mat_get(NUMBER out, mat m, int i, int j);
void mat_set(mat m, int i, int j, NUMBER x);
NUMBER *mat_ref(mat m, int i, int j);
void mat_zero(mat m);
void mat_identity(mat m);
void mat_copy(mat to, mat from);
void mat_clear(mat m);
int mat_same(mat m1, mat m2);
static void mat_multiply_outofplace(mat_workspace *ws, mat out, mat in1, mat in2);
void mat_multiply(mat_workspace *ws, mat out, mat in1, mat in2);
void mat_det(mat_workspace *ws, NUMBER out, mat in);
static void mat_pseudoinverse_outofplace(mat_workspace *ws, mat out, mat in);
void mat_pseudoinverse(mat_workspace *ws, mat out, mat in);
void mat_trace(NUMBER out, mat in);
void mat_print(mat in);

#endif