diff --git a/qext.c b/qext.c
new file mode 100644
index 0000000..fe6933e
--- /dev/null
+++ b/qext.c
@@ -0,0 +1,131 @@
+#include "qext.h"
+
+const int qext_coefficients[QEXT_RANK+1] = {36, 0, -8, 0, 1};
+
+#include <stdio.h>
+#include <malloc.h>
+
+#define LOOP(i,n) for(int i = 0; i < (n); i++)
+
+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_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)
+{
+	static initialized = 0;
+	static mpq_t result[2*x->rk-1];
+	static mpq_t coefficient[x->rk];
+	static mpq_t tmp;
+
+	if(!initialized) {
+		run_before = 1;
+		mpq_init(tmp);
+		LOOP(i, 2*x->rk-1) {
+			mpq_init(result[i]);
+		}
+		LOOP(i, x->rk) {
+			mpq_init(coefficient[i]);
+			mpq_set_si(coefficient[i], -qext_coefficients[i], qext_coefficients[x->rk]);
+		}
+	}
+
+	// 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, 2*x->rk-1) mpq_set_ui(result[i], 0, 1);
+
+	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, 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]);
+	LOOP(i, x->rk) mpq_clear(coefficient[i]);
+	mpq_clear(tmp);
+	*/
+}
+
+void qext_div(qext_number out, qext_number x, qext_number y);
diff --git a/qext.h b/qext.h
new file mode 100644
index 0000000..2e81d32
--- /dev/null
+++ b/qext.h
@@ -0,0 +1,37 @@
+#ifndef QEXT_H
+#define QEXT_H
+
+#include <gmp.h>
+
+#define QEXT_RANK 4
+
+struct qext_number_internal {
+	int rk;
+	mpq_t *a;
+};
+
+typedef struct qext_number_internal qext_number[1];
+
+void qext_init(qext_number x);
+void qext_clear(qext_number x);
+void qext_set(qext_number x, qext_number y);
+void qext_set_int(qext_number x, int y);
+void qext_add(qext_number result, qext_number x, qext_number y);
+void qext_sub(qext_number result, qext_number x, qext_number y);
+void qext_mul(qext_number out, qext_number x, qext_number y);
+void qext_div(qext_number out, qext_number x, qext_number y);
+int qext_cmp(qext_number x, qext_number y);
+void qext_print(qext_number x);
+void qext_sprint(qext_number x);
+
+static void qext_set_zero(qext_number x)
+{
+	qext_set_int(x, 0);
+}
+
+static void qext_set_one(qext_number x)
+{
+	qext_set_int(x, 1);
+}
+
+#endif