#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)

#ifdef QEXT_SQRT5

typedef mpq_t NUMBER [2];
#define INIT(x)         do { mpq_init((x)[0]);         mpq_init((x)[1]);       } while(0)
#define CLEAR(x)        do { mpq_clear((x)[0]);        mpq_clear((x)[1]);      } while(0)
#define SET(x,y)        do { mpq_set((x)[0],(y)[0]);   mpq_set((x)[1],(y)[1]); } while(0)
#define SET_INT(x,y)    do { mpq_set_si((x)[0],(y),1); mpq_set_si((x)[1],0,1); } while(0)
#define SET_ZERO(x) SET_INT(x,0)
#define SET_ONE(x)  SET_INT(x,1)
#define ADD(x,y,z)      do { mpq_add((x)[0],(y)[0],(z)[0]); mpq_add((x)[1],(y)[1],(z)[1]); } while(0)
#define SUB(x,y,z)      do { mpq_sub((x)[0],(y)[0],(z)[0]); mpq_sub((x)[1],(y)[1],(z)[1]); } while(0)
#define MUL multiply_sqrt5
#define DIV divide_sqrt5
#define CMP compare_sqrt5
#define PRINT(x) gmp_printf("%Qd + %Qd*sqrt(5)", (x)[0], (x)[1])
#define SPRINT(buf, x) gmp_sprintf(buf, "%Qd+%Qd*sqrt(5)", (x)[0], (x)[1])

static void multiply_sqrt5(NUMBER out, NUMBER a, NUMBER b)
{
	// could be in place!!!

	mpq_t tmp, result[2];
	mpq_inits(tmp, result[0], result[1], 0);

	mpq_mul(result[0], a[1], b[1]);
	mpq_set_si(tmp, 5, 1);
	mpq_mul(result[0], result[0], tmp);
	mpq_mul(tmp, a[0], b[0]);
	mpq_add(result[0], result[0], tmp);

	mpq_mul(result[1], a[1], b[0]);
	mpq_mul(tmp, a[0], b[1]);
	mpq_add(result[1], result[1], tmp);

	mpq_set(out[0], result[0]);
	mpq_set(out[1], result[1]);

	mpq_clears(tmp, result[0], result[1], 0);
}

static void divide_sqrt5(NUMBER out, NUMBER a, NUMBER b)
{
	mpq_t denom, num, tmp, neg5, result[2];
	mpq_inits(denom, num, tmp, neg5, result[0], result[1], 0);
	mpq_set_si(neg5, -5, 1);

	mpq_mul(denom, b[1], b[1]);
	mpq_mul(denom, denom, neg5);
	mpq_mul(tmp, b[0], b[0]);
	mpq_add(denom, denom, tmp);

	mpq_mul(num, a[1], b[1]);
	mpq_mul(num, num, neg5);
	mpq_mul(tmp, a[0], b[0]);
	mpq_add(num, num, tmp);

	mpq_div(result[0], num, denom);

	mpq_mul(num, a[1], b[0]);
	mpq_mul(tmp, a[0], b[1]);
	mpq_sub(num, num, tmp);

	mpq_div(result[1], num, denom);

	mpq_set(out[0], result[0]);
	mpq_set(out[1], result[1]);

	mpq_clears(denom, num, tmp, neg5, result[0], result[1], 0);
}

static int compare_sqrt5(NUMBER x, NUMBER y)
{
	int result;
	mpq_t p, q, p2, q2, c5;
	mpq_inits(p, q, p2, q2, c5, 0);

	mpq_sub(p, x[0], y[0]);
	mpq_sub(q, x[1], y[1]);

	/*
	  want to know if p + sqrt(5) q > 0
	  if p>0 and q>0: always true
	  if p>0 and q<0: equivalent to |p|^2 > 5 |q|^2
	  if p<0 and q>0: equivalent to |p|^2 < 5 |q|^2
	  if p<0 and q<0: always false
	*/

	if(mpq_sgn(p) > 0 && mpq_sgn(q) > 0) {
		result = 1;
		goto done;
	}

	if(mpq_sgn(p) < 0 && mpq_sgn(q) < 0) {
		result = -1;
		goto done;
	}

	mpq_mul(p2, p, p);
	mpq_mul(q2, q, q);
	mpq_set_si(c5, 5, 1);
	mpq_mul(q2, q2, c5);

	if(mpq_sgn(p) > 0)
		result = mpq_cmp(p2, q2);  // this can't be zero, or else |p/q| = sqrt(5), but it's irrational
	else if(mpq_sgn(p) < 0)
		result = -mpq_cmp(p2, q2); // this can be zero if p = 0 and q = 0
	else // p = 0
		return mpq_sgn(q);

done:
	mpq_clears(p, q, p2, q2, c5, 0);
	return result;
}

#endif

#ifdef QEXT_TRIVIAL

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

#endif

#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