commit b144b7720670e644e23afd5cdebcdae5fdbd72c6 Author: Florian Stecker Date: Tue Dec 5 23:06:18 2017 +0100 triangle group 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