triangle_reflection_complex/qext.c

#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++)
#define RANK(x) ((x)->type->rank)

const int qext_type_trivial_coeffs[] = {-1, 1};
struct qext_type qext_type_trivial_v = {1, qext_type_trivial_coeffs, NULL};
struct qext_type *QT_TRIVIAL = &qext_type_trivial_v;

const int qext_type_sqrt5_coeffs[] = {-5, 0, 1};
struct qext_type qext_type_sqrt5_v = {2, qext_type_sqrt5_coeffs, NULL};
struct qext_type *QT_SQRT5 = &qext_type_sqrt5_v;

static void qext_init_type(struct qext_type *type)
{
	type->coeffs = malloc(type->rank*sizeof(mpq_t));
	LOOP(i, type->rank) {
		mpq_init(type->coeffs[i]);
		mpq_set_si(type->coeffs[i],
		           -type->integer_coeffs[i],
		           type->integer_coeffs[type->rank]);
	}
}

struct qext_type *qext_newtype(int rank, const int *coeffs)
{
	struct qext_type *type = malloc(sizeof(struct qext_type));
	type->rank = rank;
	type->integer_coeffs = coeffs;
	qext_init_type(type);

	return type;
}

void qext_init(qext_number x, struct qext_type *type)
{
	if(type->coeffs == NULL) // uninitialized default type
		qext_init_type(type);

	x->type = type;
	x->a = malloc(RANK(x)*sizeof(mpq_t));
	LOOP(i, RANK(x)) mpq_init(x->a[i]);
}

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

void qext_set(qext_number x, qext_number y)
{
	LOOP(i, RANK(x)) 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 < RANK(x); 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 < RANK(x); i++)
		mpq_set_ui(x->a[i], 0, 1);
}

void qext_add(qext_number result, qext_number x, qext_number y)
{
	LOOP(i, RANK(x)) 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, RANK(x)) mpq_sub(result->a[i], x->a[i], y->a[i]);
}

void qext_neg(qext_number result, qext_number x)
{
	LOOP(i, RANK(x)) mpq_neg(result->a[i], x->a[i]);
}

void qext_print(qext_number x)
{
	LOOP(i, RANK(x)) {
		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, RANK(x)) {
		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, RANK(x)) {
		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*RANK(x)-1];
	mpq_t tmp;

	mpq_init(tmp);
	LOOP(i, 2*RANK(x)-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, RANK(x)) LOOP(j, RANK(x)) {
		mpq_mul(tmp, x->a[i], y->a[j]);
		mpq_add(result[i+j], result[i+j], tmp);
	}

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

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

	LOOP(i, 2*RANK(x)-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
}