#include "qext.h"

int qext_rank;
mpq_t *qext_coefficient;

#include <stdio.h>
#include <malloc.h>

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

void qext_setup(int rank, const int *coeffs)
{
	qext_rank = rank;
	qext_coefficient = malloc(rank*sizeof(mpq_t));
	LOOP(i, rank) {
		mpq_init(qext_coefficient[i]);
		mpq_set_si(qext_coefficient[i], -coeffs[i], coeffs[rank]);
	}
}

void qext_init(qext_number x)
{
	x->rk = qext_rank;
	x->a = malloc(x->rk*sizeof(mpq_t));
	LOOP(i, x->rk) mpq_init(x->a[i]);
}

void qext_clear(qext_number x)
{
	LOOP(i, x->rk) mpq_clear(x->a[i]);
	free(x->a);
}

void qext_set(qext_number x, qext_number y)
{
	LOOP(i, x->rk) mpq_set(x->a[i], y->a[i]);
}

void qext_set_int(qext_number x, int y)
{
	mpq_set_si(x->a[0], y, 1);
	for(int i = 1; i < x->rk; i++)
		mpq_set_ui(x->a[i], 0, 1);
}

void qext_set_q(qext_number x, mpq_t y)
{
	mpq_set(x->a[0], y);
	for(int i = 1; i < x->rk; i++)
		mpq_set_ui(x->a[i], 0, 1);
}

void qext_add(qext_number result, qext_number x, qext_number y)
{
	LOOP(i, x->rk) mpq_add(result->a[i], x->a[i], y->a[i]);
}

void qext_sub(qext_number result, qext_number x, qext_number y)
{
	LOOP(i, x->rk) mpq_sub(result->a[i], x->a[i], y->a[i]);
}

void qext_print(qext_number x)
{
	LOOP(i, x->rk) {
		if(i == 0)
			gmp_printf("%Qd", x->a[i]);
		else if(i == 1)
			gmp_printf(" + %Qd*a", x->a[i]);
		else
			gmp_printf(" + %Qd*a^%d", x->a[i], i);
	}
}

void qext_snprint(char *str, size_t size, qext_number x)
{
	int len = 0;

	LOOP(i, x->rk) {
		if(i == 0)
			len += gmp_snprintf(str+len, size-len, "%Qd", x->a[i]);
		else if(i == 1)
			len += gmp_snprintf(str+len, size-len, " + %Qd*a", x->a[i]);
		else
			len += gmp_snprintf(str+len, size-len, " + %Qd*a^%d", x->a[i], i);
	}
}

int qext_cmp(qext_number x, qext_number y)
{
	int cmp;

	LOOP(i, x->rk) {
		cmp = mpq_cmp(x->a[i], y->a[i]);
		if(cmp)
			return cmp;
	}

	return 0;
}

void qext_mul(qext_number out, qext_number x, qext_number y)
{
	mpq_t result[2*x->rk-1];
	mpq_t tmp;

	mpq_init(tmp);
	LOOP(i, 2*x->rk-1) {
		mpq_init(result[i]);
		mpq_set_ui(result[i], 0, 1);
	}

	// rk*rk multiplications + (rk-1)*rk multiplications with fixed coefficients
	// degree 1: 1+0
	// degree 2: 4+2
	// degree 4: 16+12, including 6 times 0*

	LOOP(i, x->rk) LOOP(j, x->rk) {
		mpq_mul(tmp, x->a[i], y->a[j]);
		mpq_add(result[i+j], result[i+j], tmp);
	}

	for(int i = x->rk-2; i >= 0; i--) {
		LOOP(j, x->rk) {
			mpq_mul(tmp, qext_coefficient[j], result[x->rk+i]);
			mpq_add(result[i+j], result[i+j], tmp);
		}
	}

	LOOP(i, x->rk) mpq_set(out->a[i], result[i]);

	LOOP(i, 2*x->rk-1) mpq_clear(result[i]);
	mpq_clear(tmp);
}

void qext_div(qext_number out, qext_number x, qext_number y)
{
	// todo: implement this by solving a linear system
}