triangle group

Makefile

@@ -0,0 +1,25 @@
+HEADERS=triangle.h linalg.h queue.h
+
+#SPECIAL_OPTIONS=-O0 -g -D_DEBUG
+#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
+SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
+#SPECIAL_OPTIONS=
+
+OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
+
+all: triangle_group
+
+triangle_group: triangle.o linalg.o main.o
+	gcc $(OPTIONS) -o triangle_group triangle.o linalg.o main.o -lm -lgsl -lcblas
+
+main.o: main.c $(HEADERS)
+	gcc $(OPTIONS) -c main.c
+
+linalg.o: linalg.c $(HEADERS)
+	gcc $(OPTIONS) -c linalg.c
+
+triangle.o: triangle.c $(HEADERS)
+	gcc $(OPTIONS) -c triangle.c
+
+clean:
+	rm -f triangle_group triangle.o linalg.o main.o

linalg.c

@@ -0,0 +1,138 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <memory.h>
+
+#include "linalg.h"
+
+#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
+#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
+
+/*********************************************** temporary storage ********************************************************/
+
+workspace_t *workspace_alloc(int n)
+{
+  workspace_t *result = (workspace_t*)malloc(sizeof(workspace_t));
+  result->n = n;
+  result->work_nonsymmv = gsl_eigen_nonsymmv_alloc(n);
+  result->work_sv = gsl_vector_alloc(n);
+  result->eval_complex = gsl_vector_complex_alloc(n);
+  result->evec_complex = gsl_matrix_complex_alloc(n, n);
+  result->eval_real = gsl_vector_alloc(n);
+  result->evec_real = gsl_matrix_alloc(n, n);
+  result->tmp = gsl_matrix_alloc(n, n);
+  result->permutation = gsl_permutation_alloc(n);
+  for(int i = 0; i < 20; i++)
+    result->stack[i] = gsl_matrix_alloc(n, n);
+
+  return result;
+}
+
+void workspace_free(workspace_t *workspace)
+{
+  gsl_eigen_nonsymmv_free(workspace->work_nonsymmv);
+  gsl_vector_free(workspace->work_sv);
+  gsl_vector_complex_free(workspace->eval_complex);
+  gsl_matrix_complex_free(workspace->evec_complex);
+  gsl_vector_free(workspace->eval_real);
+  gsl_matrix_free(workspace->evec_real);
+  gsl_matrix_free(workspace->tmp);
+  gsl_permutation_free(workspace->permutation);
+  for(int i = 0; i < 20; i++)
+    gsl_matrix_free(workspace->stack[i]);
+}
+
+/************************************************** basic operations ********************************************************/
+
+void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws)
+{
+  int s;
+  gsl_matrix_memcpy(ws->tmp, in);
+  gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s);
+  gsl_linalg_LU_invert(ws->tmp, ws->permutation, out);
+}
+
+void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws)
+{
+  invert(conjugator, out, ws);   // use out to temporarily store inverse conjugator
+  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, ws->tmp);  // in * conjugator^{-1}
+  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, ws->tmp, 0.0, out);
+}
+
+void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out)
+{
+  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, out);
+}
+
+void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
+{
+  gsl_matrix_memcpy(ws->tmp, g);
+  gsl_linalg_SV_decomp(ws->tmp, ws->evec_real, ws->eval_real, ws->work_sv);
+
+  for(int i = 0; i < ws->n - 1; i++)
+    mu[i] = log(gsl_vector_get(ws->eval_real, i) / gsl_vector_get(ws->eval_real, i+1));
+}
+
+void initialize(gsl_matrix *g, double *data, int x, int y)
+{
+  gsl_matrix_view view = gsl_matrix_view_array(data, x, y);
+  gsl_matrix_memcpy(g, &view.matrix);
+}
+
+void rotation_matrix(gsl_matrix *g, double *vector)
+{
+  double normalized[3];
+  double norm = sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2]);
+  for(int i = 0; i < 3; i++)
+    normalized[i] = vector[i] / norm;
+  gsl_matrix_set_identity(g);
+  gsl_matrix_set(g, 0, 0, cos(norm));
+  gsl_matrix_set(g, 0, 1, -sin(norm) * normalized[2]);
+  gsl_matrix_set(g, 0, 2, +sin(norm) * normalized[1]);
+  gsl_matrix_set(g, 1, 0, +sin(norm) * normalized[2]);
+  gsl_matrix_set(g, 1, 1, cos(norm));
+  gsl_matrix_set(g, 1, 2, -sin(norm) * normalized[0]);
+  gsl_matrix_set(g, 2, 0, -sin(norm) * normalized[1]);
+  gsl_matrix_set(g, 2, 1, +sin(norm) * normalized[0]);
+  gsl_matrix_set(g, 2, 2, cos(norm));
+  for(int i = 0; i < 3; i++)
+    for(int j = 0; j < 3; j++)
+      g->data[i * g->tda + j] += (1 - cos(norm)) * normalized[i] * normalized[j];
+}
+
+double trace(gsl_matrix *g)
+{
+  return gsl_matrix_get(g, 0, 0) + gsl_matrix_get(g, 1, 1) + gsl_matrix_get(g, 2, 2);
+}
+
+double determinant(gsl_matrix *g, workspace_t *ws)
+{
+  int s;
+  gsl_matrix_memcpy(ws->tmp, g);
+  gsl_linalg_LU_decomp(ws->tmp, ws->permutation, &s);
+  return gsl_linalg_LU_det(ws->tmp, s);
+}
+
+void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws)
+{
+  gsl_eigen_nonsymmv_params(1, ws->work_nonsymmv);
+  gsl_eigen_nonsymmv(g, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv);
+  gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC);
+
+  int real = 1;
+  for(int i = 0; i < 4; i++)
+    if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0)
+      real = 0;
+
+  if(!real) {
+    printf("We have non-real eigenvalues!\n");
+    exit(1);
+  }
+
+  for(int i = 0; i < ws->n - 1; i++) {
+    mu[i] = log(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, i+1)));
+  }
+}

linalg.h

@@ -0,0 +1,40 @@
+#ifndef LINALG_H
+#define LINALG_H
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <memory.h>
+
+#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
+#define FCMP(x, y) gsl_fcmp(x, y, 1e-10)
+
+typedef struct _workspace {
+  int n;
+  gsl_eigen_nonsymmv_workspace *work_nonsymmv;
+  gsl_vector *work_sv;
+  gsl_vector_complex *eval_complex;
+  gsl_matrix_complex *evec_complex;
+  gsl_vector *eval_real;
+  gsl_matrix *evec_real;
+  gsl_matrix *tmp;
+  gsl_permutation *permutation;
+  gsl_matrix *stack[20];
+} workspace_t;
+
+workspace_t *workspace_alloc(int n);
+void workspace_free(workspace_t *workspace);
+void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws);
+void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws);
+void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out);
+void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
+void initialize(gsl_matrix *g, double *data, int x, int y);
+void rotation_matrix(gsl_matrix *g, double *vector);
+void jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws);
+double trace(gsl_matrix *g);
+double determinant(gsl_matrix *g, workspace_t *ws);
+
+#endif

main.c

@@ -0,0 +1,120 @@
+#include "triangle.h"
+#include "linalg.h"
+
+#define MAX_ELEMENTS 20000
+
+void initialize_triangle_generators(gsl_matrix **gen, double a1, double a2, double a3, double s)
+{
+  for(int i = 0; i < 3; i++)
+    gsl_matrix_set_identity(gen[i]);
+
+  gsl_matrix_set(gen[0], 0, 0, -1);
+  gsl_matrix_set(gen[0], 1, 0, 2*s*cos(M_PI*a1));
+  gsl_matrix_set(gen[0], 2, 0, 2/s*cos(M_PI*a2));
+
+  gsl_matrix_set(gen[1], 0, 1, 2/s*cos(M_PI*a1));
+  gsl_matrix_set(gen[1], 1, 1, -1);
+  gsl_matrix_set(gen[1], 2, 1, 2*s*cos(M_PI*a3));
+
+  gsl_matrix_set(gen[2], 0, 2, 2*s*cos(M_PI*a2));
+  gsl_matrix_set(gen[2], 1, 2, 2/s*cos(M_PI*a3));
+  gsl_matrix_set(gen[2], 2, 2, -1);
+}
+
+char *print_word(groupelement_t *g, char *str)
+{
+  int i = g->length - 1;
+
+  str[g->length] = 0;
+  while(g->parent) {
+    str[i--] = 'a' + g->letter;
+    g = g->parent;
+  }
+
+  return str;
+}
+
+int main(int argc, char *argv[])
+{
+  groupelement_t *group;
+  workspace_t *ws;
+  gsl_matrix **matrices;
+  gsl_matrix *gen[3];
+  gsl_matrix *tmp, *tmp2, *coxeter;
+  double mu[2];
+  char buf[100];
+
+  if(argc < 2) {
+    fprintf(stderr, "Too few arguments!\n");
+    return 1;
+  }
+
+  // allocate stuff
+  group = malloc(MAX_ELEMENTS*sizeof(groupelement_t));
+  ws = workspace_alloc(3);
+  matrices = malloc(MAX_ELEMENTS*sizeof(gsl_matrix*));
+  for(int i = 0; i < MAX_ELEMENTS; i++)
+    matrices[i] = gsl_matrix_alloc(3, 3);
+  for(int i = 0; i < 3; i++)
+    gen[i] = gsl_matrix_alloc(3, 3);
+  tmp = gsl_matrix_alloc(3, 3);
+  tmp2 = gsl_matrix_alloc(3, 3);
+  coxeter = gsl_matrix_alloc(3, 3);
+
+  // the real action
+  generate_triangle_group(group, MAX_ELEMENTS, 5, 9, 7);
+  //  initialize_triangle_generators(gen, 1.0/5, 1.0/5, 1.0/5, 1.0); // Fuchsian
+  initialize_triangle_generators(gen, 3.0/5.0, 5.0/9.0, 3.0/7.0, atof(argv[1]));
+
+  gsl_matrix_set_identity(matrices[0]);
+  for(int i = 1; i < MAX_ELEMENTS; i++)
+    multiply(matrices[group[i].parent->id], gen[group[i].letter], matrices[i]);
+
+  for(int i = 0; i < MAX_ELEMENTS; i++) {
+    cartan_calc(matrices[i], mu, ws);
+    printf("%d %f %f 0x0000ff %s\n", group[i].length, mu[0], mu[1], print_word(&group[i], buf));
+    //invert(matrices[i], tmp, ws);
+    //printf("%d %f %f %f 0x0000ff\n", group[i].length, determinant(matrices[i], ws)*trace(matrices[i]), determinant(matrices[i], ws)*trace(tmp), determinant(matrices[i], ws));
+  }
+
+  /*
+  for(int p = 0; p < 6; p++) {
+    groupelement_t *cur = &group[0];
+    while(1) {
+      if(!(cur = cur->adj[p%3]))
+	break;
+
+      cartan_calc(matrices[cur->id], mu, ws);
+      printf("%d %f %f 0xff0000\n", group[cur->id].length, mu[0], mu[1]);
+
+      if(!(cur = cur->adj[(p+p/3+1)%3]))
+	break;
+
+      cartan_calc(matrices[cur->id], mu, ws);
+      printf("%d %f %f 0xff0000\n", group[cur->id].length, mu[0], mu[1]);
+
+      if(!(cur = cur->adj[(p-p/3+2)%3]))
+	break;
+
+      cartan_calc(matrices[cur->id], mu, ws);
+      printf("%d %f %f 0xff0000\n", group[cur->id].length, mu[0], mu[1]);
+      //      invert(matrices[cur->id], tmp, ws);
+      //      printf("%d %f %f %f 0xff0000\n", group[cur->id].length, determinant(matrices[cur->id], ws)*trace(matrices[cur->id]), determinant(matrices[cur->id], ws)*trace(tmp), determinant(matrices[cur->id], ws));
+    }
+  }
+  */
+
+  // free stuff
+  for(int i = 0; i < MAX_ELEMENTS; i++)
+    gsl_matrix_free(matrices[i]);
+  for(int i = 0; i < 3; i++)
+    gsl_matrix_free(gen[i]);
+  gsl_matrix_free(tmp);
+  gsl_matrix_free(tmp2);
+  gsl_matrix_free(coxeter);
+  free(matrices);
+  free(group);
+  workspace_free(ws);
+
+  return 0;
+}

queue.h

@@ -0,0 +1,46 @@
+#ifndef QUEUE_H
+#define QUEUE_H
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#define QUEUE_SIZE 50000
+
+#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
+#ifdef _DEBUG
+#define LOG(msg, ...) fprintf(stderr, msg, ##__VA_ARGS__)
+#else
+#define LOG(msg, ...)
+#endif
+
+typedef struct {
+  unsigned int start;
+  unsigned int end;
+  int data[QUEUE_SIZE];
+} queue_t;
+
+static void queue_init(queue_t *q)
+{
+  q->start = q->end = 0;
+}
+
+static void queue_put(queue_t *q, int x)
+{
+  q->data[q->end++] = x;
+  q->end %= QUEUE_SIZE;
+
+  ERROR(q->start == q->end, "The queue is full! Increase QUEUE_SIZE\n");
+}
+
+static int queue_get(queue_t *q)
+{
+  if(q->start == q->end)
+    return -1;
+
+  int result = q->data[q->start++];
+  q->start %= QUEUE_SIZE;
+
+  return result;
+}
+
+#endif

triangle.c

@@ -0,0 +1,61 @@
+#include <stdio.h>
+#include <string.h>
+#include <memory.h>
+
+#include "triangle.h"
+#include "queue.h"
+
+static groupelement_t *find_adj(groupelement_t *cur, int gen, int other, int order)
+{
+	groupelement_t *next = cur->adj[other];
+	for(int i = 0; i < order - 1; i++) {
+		if(!next)
+			break;
+		next = next->adj[gen];
+		if(!next)
+			break;
+		next = next->adj[other];
+	}
+	return next;
+}
+
+int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3)
+{
+	queue_t q;
+	int id, n;
+	groupelement_t *cur;
+
+	int k[3] = {k1, k2, k3};
+
+	memset(group, 0, size*sizeof(groupelement_t));
+	n = 1;
+	queue_init(&q);
+	queue_put(&q, 0);
+
+	while((id = queue_get(&q)) != -1) {
+		cur = &group[id];
+
+		for(int gen = 0; gen < 3; gen++) {
+			if(!cur->adj[gen])
+				cur->adj[gen] = find_adj(cur, gen, (gen+2)%3, k[(gen+1)%3]);
+
+			if(!cur->adj[gen])
+				cur->adj[gen] = find_adj(cur, gen, (gen+1)%3, k[(gen+2)%3]);
+
+			if(!cur->adj[gen] && n < size) {
+				cur->adj[gen] = &group[n];
+				cur->adj[gen]->id = n;
+				cur->adj[gen]->parent = cur;
+				cur->adj[gen]->letter = gen;
+				cur->adj[gen]->length = cur->length + 1;
+				queue_put(&q, n);
+				n++;
+			}
+
+			if(cur->adj[gen])
+				cur->adj[gen]->adj[gen] = cur;
+		}
+	}
+
+	return n;
+}

triangle.h

@@ -0,0 +1,20 @@
+#ifndef TRIANGLE_H
+#define TRIANGLE_H
+
+#include <stdio.h>
+#include <string.h>
+#include <memory.h>
+
+#include "queue.h"
+
+typedef struct _groupelement {
+  int id;
+  int length;
+  struct _groupelement *adj[3];
+  struct _groupelement *parent;
+  int letter;
+} groupelement_t;
+
+int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);
+
+#endif

triangle_group.plt

@@ -0,0 +1,20 @@
+#if(!exists("param")) param = log(0.3810364354550512)
+if(!exists("param")) param = log(2.624420939708161)
+
+print sprintf("param = %f", exp(param))
+file = sprintf("< ./triangle_group %f", exp(param))
+
+set zeroaxis
+set size square
+set xrange [0:20]
+set yrange [0:20]
+set grid
+set parametric
+
+plot file using 1:2 w p pt 7 ps 1 lc 1 t ""
+
+pause mouse keypress
+if(MOUSE_KEY == 60) param=param-0.02
+if(MOUSE_KEY == 62) param=param+0.02
+if(MOUSE_KEY == 13) param=0
+if(MOUSE_KEY != 113) reread