working version of Q extensions

Makefile

@@ -1,8 +1,8 @@
-HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h
+HEADERS=linalg.h mat.h coxeter.h enumerate_triangle_group.h parallel.h qext.h
 
-#SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT_SQRT5
+SPECIAL_OPTIONS=-O0 -g -D_DEBUG -DQEXT
 #SPECIAL_OPTIONS=-O3 -pg -g -funroll-loops -fno-inline
-SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT_SQRT5
+#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -DQEXT
 #SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline -mavx512f -mavx512cd -mavx512er -mavx512pf # KNL
 #SPECIAL_OPTIONS=
 
@@ -19,8 +19,8 @@ billiard_words: billiard_words.hs
 singular_values: singular_values.o coxeter.o mat.o enumerate_triangle_group.o parallel.o
 	mpicc $(OPTIONS) -o singular_values coxeter.o singular_values.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps
 
-singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o
+singular_values_barbot: singular_values_barbot.o coxeter.o mat.o enumerate_triangle_group.o parallel.o qext.o
-	mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o -lm -lgmp -lmps
+	mpicc $(OPTIONS) -o singular_values_barbot coxeter.o singular_values_barbot.o mat.o enumerate_triangle_group.o parallel.o qext.o -lm -lgmp -lmps
 
 #singular_values_mpi: singular_values_mpi.o coxeter.o mat.o
 #	mpicc $(OPTIONS) -o singular_values_mpi coxeter.o singular_values_mpi.o mat.o -lm -lgmp -lmps
@@ -55,5 +55,8 @@ mat.o: mat.c $(HEADERS)
 parallel.o: parallel.c $(HEADERS)
 	gcc $(OPTIONS) -c parallel.c
 
+qext.o: qext.c $(HEADERS)
+	gcc $(OPTIONS) -c qext.c
+
 clean:
-	rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.o singular_values_mpi.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words enumerate_triangle_group.o parallel.o
+	rm -f singular_values singular_values_barbot special_element singular_values_mpi coxeter.o linalg.o singular_values.o singular_values_barbot.o singular_values_mpi.o mat.o special_element.o convert.hi convert.o convert billiard_words.hi billiard_words.o billiard_words enumerate_triangle_group.o parallel.o qext.o

enumerate_triangle_group.c

@@ -129,33 +129,38 @@ void discriminant(NUMBER out, NUMBER x, NUMBER y)
 	CLEAR(x2);CLEAR(x3);CLEAR(y2);CLEAR(y3);CLEAR(tmp);
 }
 
-void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER s, NUMBER q)
+void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, mpq_t s, mpq_t q)
 {
 	mat_workspace *ws;
 	mat r1,r2,r3;
 	NUMBER b1,c1,a2,c2,a3,b3;
-	NUMBER sinv;
+	mpq_t sinv, qinv;
 
 	ws = mat_workspace_init(3);
-	INIT(sinv);INIT(b1);INIT(c1);INIT(a2);INIT(c2);INIT(a3);INIT(b3);
+	INIT(b1);INIT(c1);INIT(a2);INIT(c2);INIT(a3);INIT(b3);
+	mpq_init(sinv);
+	mpq_init(qinv);
 	mat_init(r1, 3);
 	mat_init(r2, 3);
 	mat_init(r3, 3);
 
 	// sinv = s^{-1}
-	SET_INT(sinv, 1);
+	mpq_inv(sinv, s);
-	DIV(sinv, sinv, s);
+	mpq_inv(qinv, q);
 
 	// c1 = rho2 q, a2 = rho3 q, b3 = rho1 q, b1 = c2 = a3 = 1/q
-	MUL(c1, rho2, q);
+	SET_Q(c1, q);
-	MUL(a2, rho3, q);
+	SET_Q(a2, q);
-	MUL(b3, rho1, q);
+	SET_Q(b3, q);
+	MUL(c1, c1, rho2);
+	MUL(a2, a2, rho3);
+	MUL(b3, b3, rho1);
 	SET_INT(b1, 1);
 	SET_INT(c2, 1);
 	SET_INT(a3, 1);
-	DIV(b1, b1, q);
+	SET_Q(b1, qinv);
-	DIV(c2, c2, q);
+	SET_Q(c2, qinv);
-	DIV(a3, a3, q);
+	SET_Q(a3, qinv);
 
 	mat_zero(r1);
 	mat_zero(r2);
@@ -183,17 +188,17 @@ void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NU
 
 	// gen[0] = diag(1,s^{-1},s)
 	SET_INT(*mat_ref(gen[0], 0, 0), 1);
-	mat_set(gen[0], 1, 1, sinv);
+	SET_Q  (*mat_ref(gen[0], 1, 1), sinv);
-	mat_set(gen[0], 2, 2, s);
+	SET_Q  (*mat_ref(gen[0], 2, 2), s);
 
 	// gen[1] = diag(s,1,s^{-1})
-	mat_set(gen[1], 0, 0, s);
+	SET_Q  (*mat_ref(gen[1], 0, 0), s);
 	SET_INT(*mat_ref(gen[1], 1, 1), 1);
-	mat_set(gen[1], 2, 2, sinv);
+	SET_Q  (*mat_ref(gen[1], 2, 2), sinv);
 
-	// gen[3] = diag(s^{-1},s,1)
+	// gen[2] = diag(s^{-1},s,1)
-	mat_set(gen[2], 0, 0, sinv);
+	SET_Q  (*mat_ref(gen[2], 0, 0), sinv);
-	mat_set(gen[2], 1, 1, s);
+	SET_Q  (*mat_ref(gen[2], 1, 1), s);
 	SET_INT(*mat_ref(gen[2], 2, 2), 1);
 
 	// gen[0] = r2 * gen[0] * r3
@@ -214,7 +219,9 @@ void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NU
 	mat_pseudoinverse(ws, gen[5], gen[2]);
 
 	mat_workspace_clear(ws);
-	CLEAR(sinv);CLEAR(b1);CLEAR(c1);CLEAR(a2);CLEAR(c2);CLEAR(a3);CLEAR(b3);
+	CLEAR(b1);CLEAR(c1);CLEAR(a2);CLEAR(c2);CLEAR(a3);CLEAR(b3);
+	mpq_clear(sinv);
+	mpq_clear(qinv);
 	mat_clear(r1);
 	mat_clear(r2);
 	mat_clear(r3);
@@ -239,31 +246,29 @@ void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3
 }
 #endif
 
-#ifdef QEXT_SQRT5
+#ifdef QEXT
 void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s_, mpq_t q_)
 {
-	NUMBER rho, c, one, s, q;
+	NUMBER rho, c, one, s;
-	INIT(rho);INIT(c);INIT(one);INIT(s);INIT(q);
+	INIT(rho);INIT(c);INIT(one);INIT(s);
 
 	// c = - (1+sqrt(5))/2
-	mpq_set_si(c[0], -1, 2);
+	mpq_set_si(c->a[0], -1, 2);
-	mpq_set_si(c[1], -1, 2);
+	mpq_set_si(c->a[1], -1, 2);
 	SET_ONE(one);
 
-	mpq_set(s[0], s_);
+	SET_Q(s, s_);
-	mpq_set_ui(s[1], 0, 1);
-	mpq_set(q[0], q_);
-	mpq_set_ui(q[1], 0, 1);
 
+	// rho = s^2 + cs + 1
 	SET(rho, one);
 	MUL(rho, rho, s);
 	ADD(rho, rho, c);
 	MUL(rho, rho, s);
 	ADD(rho, rho, one);
 
-	generators_triangle_rotation_generic(gen, rho, rho, rho, s, q);
+	generators_triangle_rotation_generic(gen, rho, rho, rho, s_, q_);
 
-	CLEAR(rho);CLEAR(c);CLEAR(one);CLEAR(s);CLEAR(q);
+	CLEAR(rho);CLEAR(c);CLEAR(one);CLEAR(s);
 }
 #endif
 

enumerate_triangle_group.h

@@ -10,15 +10,15 @@
 int solve_characteristic_polynomial(mps_context *solv, mps_monomial_poly *poly, mpq_t tr, mpq_t trinv, double *eigenvalues);
 void continued_fraction_approximation(mpq_t out, double in, int level);
 
-// geneartors
+// generators
-void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, NUMBER s, NUMBER q);
+void generators_triangle_rotation_generic(mat *gen, NUMBER rho1, NUMBER rho2, NUMBER rho3, mpq_t s, mpq_t q);
 
 #ifdef QEXT_TRIVIAL
 void generators_triangle_rotation_2346_rational(mat *gen, int p1, int p2, int p3, mpq_t s, mpq_t q);
 #endif
 
-#ifdef QEXT_SQRT5
+#ifdef QEXT
-void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s, mpq_t q);
+void generators_triangle_rotation_555_barbot(mat *gen, mpq_t s_, mpq_t q_);
 #endif
 
 // general functions

mat.h

@@ -8,119 +8,25 @@
 
 // library for matrix computations in variable rings (based on GMP types)
 
-#ifdef QEXT_SQRT5
+#ifdef QEXT
 
-typedef mpq_t NUMBER [2];
+#include "qext.h"
-#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)
+#define NUMBER qext_number
-{
+#define INIT qext_init
-	// could be in place!!!
+#define CLEAR qext_clear
-
+#define SET qext_set
-	mpq_t tmp, result[2];
+#define SET_INT qext_set_int
-	mpq_inits(tmp, result[0], result[1], 0);
+#define SET_Q qext_set_q
-
+#define SET_ZERO(x) qext_set_int((x),0)
-	mpq_mul(result[0], a[1], b[1]);
+#define SET_ONE(x) qext_set_int((x),1)
-	mpq_set_si(tmp, 5, 1);
+#define ADD qext_add
-	mpq_mul(result[0], result[0], tmp);
+#define SUB qext_sub
-	mpq_mul(tmp, a[0], b[0]);
+#define MUL qext_mul
-	mpq_add(result[0], result[0], tmp);
+// #define DIV qext_div
-
+#define CMP qext_cmp
-	mpq_mul(result[1], a[1], b[0]);
+#define PRINT qext_print
-	mpq_mul(tmp, a[0], b[1]);
+#define SNPRINT qext_snprint
-	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
 

qext.c

@@ -1,15 +1,26 @@
 #include "qext.h"
 
-const int qext_coefficients[QEXT_RANK+1] = {36, 0, -8, 0, 1};
+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->rk = qext_rank;
 	x->a = malloc(x->rk*sizeof(mpq_t));
 	LOOP(i, x->rk) mpq_init(x->a[i]);
 }
@@ -32,6 +43,13 @@ void qext_set_int(qext_number x, int y)
 		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]);
@@ -83,21 +101,13 @@ int qext_cmp(qext_number x, qext_number y)
 
 void qext_mul(qext_number out, qext_number x, qext_number y)
 {
-	static initialized = 0;
+	mpq_t result[2*x->rk-1];
-	static mpq_t result[2*x->rk-1];
+	mpq_t tmp;
-	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]);
-		}
+		mpq_set_ui(result[i], 0, 1);
-		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
@@ -105,8 +115,6 @@ void qext_mul(qext_number out, qext_number x, qext_number y)
 	// 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);
@@ -114,18 +122,18 @@ void qext_mul(qext_number out, qext_number x, qext_number y)
 
 	for(int i = x->rk-2; i >= 0; i--) {
 		LOOP(j, x->rk) {
-			mpq_mul(tmp, coefficient[j], result[x->rk+i]);
+			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]);
-	LOOP(i, x->rk) mpq_clear(coefficient[i]);
 	mpq_clear(tmp);
-	*/
 }
 
-void qext_div(qext_number out, qext_number x, qext_number y);
+void qext_div(qext_number out, qext_number x, qext_number y)
+{
+	// todo: implement this by solving a linear system
+}

qext.h

@@ -3,8 +3,6 @@
 
 #include <gmp.h>
 
-#define QEXT_RANK 4
-
 struct qext_number_internal {
 	int rk;
 	mpq_t *a;
@@ -12,26 +10,18 @@ struct qext_number_internal {
 
 typedef struct qext_number_internal qext_number[1];
 
+void qext_setup(int rank, const int *coeffs);
 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_set_q(qext_number x, mpq_t 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);
+void qext_snprint(char *str, size_t len, 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

singular_values_barbot.c

@@ -145,12 +145,15 @@ int invariants_trace_loxodromic(group_t *group, mat *matrices, struct result **i
 
 	// check if loxodromic
 	NUMBER disc, zero;
+	double disc_real;
 	INIT(disc);
 	INIT(zero);
 	SET_ZERO(zero);
 	for(int i = 0; i < nuniq; i++) {
 		discriminant(disc, invariants[i]->tr, invariants[i]->trinv);
-		invariants[i]->disc_sign = CMP(disc, zero);
+		disc_real = mpq_get_d(disc->a[0]) + sqrt(5)*mpq_get_d(disc->a[1]);
+		gmp_printf("%Qd %Qd %f\n", disc->a[0], disc->a[1], disc_real);
+		invariants[i]->disc_sign = disc_real > 0 ? 1 : disc_real < 0 ? -1 : 0;
 	}
 	CLEAR(disc);
 	CLEAR(zero);
@@ -417,6 +420,9 @@ int main(int argc, char *argv[])
 	g->qend = atoi(argv[9]);
 	g->qdenom = atoi(argv[10]);
 
+	const int coefficients[] = {-5, 0, 1};
+	qext_setup(2, coefficients);
+
 	// initialize
 	parallel_context *ctx = parallel_init();
 	parallel_set_datasize_and_callbacks(ctx, init_node, do_computation, destroy_node,
@@ -537,13 +543,13 @@ int main(int argc, char *argv[])
 			mpq_clear(tmp);
 			*/
 
-			SPRINT(buf, n.distinct_invariants[i]->tr);
+			SNPRINT(buf, sizeof(buf), n.distinct_invariants[i]->tr);
-			SPRINT(buf2, n.distinct_invariants[i]->trinv);
+			SNPRINT(buf2, sizeof(buf2), n.distinct_invariants[i]->trinv);
-			printf("%d %s %f %f %d\n",
+			printf("%d %s %d\n",
 			       n.distinct_invariants[i]->id,
 			       print_word(&n.group->elements[n.distinct_invariants[i]->id], buf3),
-			       mpq_get_d(n.distinct_invariants[i]->tr[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->tr[1]),
+//			       mpq_get_d(n.distinct_invariants[i]->tr[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->tr[1]),
-			       mpq_get_d(n.distinct_invariants[i]->trinv[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->trinv[1]),
+//			       mpq_get_d(n.distinct_invariants[i]->trinv[0])+sqrt(5)*mpq_get_d(n.distinct_invariants[i]->trinv[1]),
 			       n.distinct_invariants[i]->disc_sign);
 
 		}