compute complex traces

qext.c

@@ -1,5 +1,7 @@
 #include "qext.h"
 
+#include <assert.h>
+
 int qext_rank;
 mpq_t *qext_coefficient;
 
@@ -9,13 +11,20 @@ mpq_t *qext_coefficient;
 #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;
+struct qext_type *QT_TRIVIAL = &(struct qext_type) {
+	.rank = 1,
+	.integer_coeffs = (const int[]) {-1, 1}
+};
 
-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;
+struct qext_type *QT_SQRT5 = &(struct qext_type) {
+	.rank = 2,
+	.integer_coeffs = (const int[]) {-5, 0, 1}
+};
+
+struct qext_type *QT_GAUSS_SQRT5 = &(struct qext_type) {
+	.rank = 4,
+	.integer_coeffs = (const int[]) {36, 0, -8, 0, 1}
+};
 
 static void qext_init_type(struct qext_type *type)
 {
@@ -161,7 +170,86 @@ void qext_mul(qext_number out, qext_number x, qext_number y)
 	mpq_clear(tmp);
 }
 
+void qext_inv(qext_number out, qext_number in)
+{
+	qext_number one;
+	qext_init(one, in->type);
+	qext_set_int(one, 1);
+	qext_div(out, one, in);
+	qext_clear(one);
+}
+
 void qext_div(qext_number out, qext_number x, qext_number y)
 {
-	// todo: implement this by solving a linear system
+	int rk = RANK(x);
+	struct qext_type *type = x->type;
+
+	mpq_t tmp;
+	mpq_t result[rk];
+	mpq_t matrix[rk*rk];
+
+	mpq_init(tmp);
+	LOOP(i, rk) mpq_init(result[i]);
+	LOOP(i, rk*rk) mpq_init(matrix[i]);
+
+	// initialize a matrix which represents multiplication by y
+	// that means the i-th column is y*a^i
+	LOOP(i, rk)
+		mpq_set(matrix[i*rk], y->a[i]);
+	LOOP(j, rk-1) { // columns
+		mpq_mul(matrix[j+1], matrix[(rk-1)*rk+j], type->coeffs[0]);
+		LOOP(i, rk-1) { // rows
+			mpq_mul(matrix[(i+1)*rk+(j+1)], matrix[(rk-1)*rk+j], type->coeffs[i+1]);
+			mpq_add(matrix[(i+1)*rk+(j+1)], matrix[(i+1)*rk+(j+1)], matrix[i*rk+j]);
+		}
+	}
+
+	// initialize result as x
+	LOOP(i, rk) mpq_set(result[i], x->a[i]);
+
+	// use Gaussian elimination to solve the system matrix * out = result
+	LOOP(j, rk)
+	{
+		// find nonzero entry
+		int row = j;
+		while(row < rk && mpq_sgn(matrix[row*rk+j]) == 0) row++;
+
+		// if the entire column is zero, the matrix was not invertible
+		// this could happen if the polynomial is not irreducible / we're not in a field
+		assert(row != rk);
+
+		// permute rows
+		if(row != j) {
+			mpq_swap(result[row], result[j]);
+			LOOP(k, rk) {
+				mpq_swap(matrix[row*rk+k], matrix[j*rk+k]);
+			}
+		}
+
+		// normalize
+		LOOP(k, rk)
+			if(k != j)
+				mpq_div(matrix[j*rk+k], matrix[j*rk+k], matrix[j*rk+j]);
+		mpq_div(result[j], result[j], matrix[j*rk+j]);
+		mpq_set_ui(matrix[j*rk+j], 1, 1);
+
+		// subtract
+		LOOP(i, rk) {
+			if(i == j)
+				continue;
+			mpq_mul(tmp, matrix[i*rk+j], result[j]);
+			mpq_sub(result[i], result[i], tmp);
+			for(int k = j+1; k < rk; k++) {
+				mpq_mul(tmp, matrix[i*rk+j], matrix[j*rk+k]);
+				mpq_sub(matrix[i*rk+k], matrix[i*rk+k], tmp);
+			}
+			mpq_set_ui(matrix[i*rk+j], 0, 1); // it isn't strictly necessary to do this as we won't use this entry again, but helps debugging
+		}
+	}
+
+	LOOP(i, rk) mpq_set(out->a[i], result[i]);
+
+	mpq_clear(tmp);
+	LOOP(i, rk) mpq_clear(result[i]);
+	LOOP(i, rk*rk) mpq_clear(matrix[i]);
 }