From b144b7720670e644e23afd5cdebcdae5fdbd72c6 Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Tue, 5 Dec 2017 23:06:18 +0100 Subject: [PATCH] triangle group --- Makefile | 25 ++++++++ linalg.c | 138 +++++++++++++++++++++++++++++++++++++++++++++ linalg.h | 40 +++++++++++++ main.c | 120 +++++++++++++++++++++++++++++++++++++++ queue.h | 46 +++++++++++++++ triangle.c | 61 ++++++++++++++++++++ triangle.h | 20 +++++++ triangle_group.plt | 20 +++++++ 8 files changed, 470 insertions(+) create mode 100644 Makefile create mode 100644 linalg.c create mode 100644 linalg.h create mode 100644 main.c create mode 100644 queue.h create mode 100644 triangle.c create mode 100644 triangle.h create mode 100644 triangle_group.plt diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..7e0e4b5 --- /dev/null +++ b/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 diff --git a/linalg.c b/linalg.c new file mode 100644 index 0000000..2514699 --- /dev/null +++ b/linalg.c @@ -0,0 +1,138 @@ +#include +#include +#include +#include +#include +#include +#include + +#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))); + } +} diff --git a/linalg.h b/linalg.h new file mode 100644 index 0000000..aac7df9 --- /dev/null +++ b/linalg.h @@ -0,0 +1,40 @@ +#ifndef LINALG_H +#define LINALG_H + +#include +#include +#include +#include +#include +#include +#include + +#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 diff --git a/main.c b/main.c new file mode 100644 index 0000000..3278513 --- /dev/null +++ b/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; +} diff --git a/queue.h b/queue.h new file mode 100644 index 0000000..2be1af0 --- /dev/null +++ b/queue.h @@ -0,0 +1,46 @@ +#ifndef QUEUE_H +#define QUEUE_H + +#include +#include + +#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 diff --git a/triangle.c b/triangle.c new file mode 100644 index 0000000..011d778 --- /dev/null +++ b/triangle.c @@ -0,0 +1,61 @@ +#include +#include +#include + +#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; +} diff --git a/triangle.h b/triangle.h new file mode 100644 index 0000000..d78241f --- /dev/null +++ b/triangle.h @@ -0,0 +1,20 @@ +#ifndef TRIANGLE_H +#define TRIANGLE_H + +#include +#include +#include + +#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 diff --git a/triangle_group.plt b/triangle_group.plt new file mode 100644 index 0000000..0f90f46 --- /dev/null +++ b/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