From 0c07494110f7c7390d97f3ba15bebd0f1c6d501b Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 1 Jul 2020 09:59:16 -0500 Subject: [PATCH] compute jordan projection for triangle rotation group --- Makefile | 25 ++++ linalg.c | 353 ++++++++++++++++++++++++++++++++++++++++++++ linalg.h | 98 ++++++++++++ queue.h | 46 ++++++ singular_values.c | 210 ++++++++++++++++++++++++++ singular_values.plt | 23 +++ triangle.c | 72 +++++++++ triangle.h | 21 +++ 8 files changed, 848 insertions(+) create mode 100644 Makefile create mode 100644 linalg.c create mode 100644 linalg.h create mode 100644 queue.h create mode 100644 singular_values.c create mode 100644 singular_values.plt create mode 100644 triangle.c create mode 100644 triangle.h diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..fa862ff --- /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: singular_values + +singular_values: singular_values.o triangle.o linalg.o + gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o -lm -lgsl -lcblas + +singular_values.o: singular_values.c $(HEADERS) + gcc $(OPTIONS) -c singular_values.c + +linalg.o: linalg.c $(HEADERS) + gcc $(OPTIONS) -c linalg.c + +triangle.o: triangle.c $(HEADERS) + gcc $(OPTIONS) -c triangle.c + +clean: + rm -f singular_values triangle.o linalg.o singular_values.o diff --git a/linalg.c b/linalg.c new file mode 100644 index 0000000..680c6a6 --- /dev/null +++ b/linalg.c @@ -0,0 +1,353 @@ +#include +#include +#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_symmv = gsl_eigen_symmv_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->permutation = gsl_permutation_alloc(n); + + result->tmp_mat = malloc(MAX_TEMP_MATRICES*sizeof(gsl_matrix*)); + for(int i = 0; i < MAX_TEMP_MATRICES; i++) + result->tmp_mat[i] = gsl_matrix_alloc(3, 3); + result->tmp_mat_used = 0; + result->tmp_vec = malloc(MAX_TEMP_MATRICES*sizeof(gsl_vector*)); + for(int i = 0; i < MAX_TEMP_MATRICES; i++) + result->tmp_vec[i] = gsl_vector_alloc(3); + result->tmp_vec_used = 0; + return result; +} + +void workspace_free(workspace_t *workspace) +{ + gsl_eigen_nonsymmv_free(workspace->work_nonsymmv); + gsl_eigen_symmv_free(workspace->work_symmv); + 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_permutation_free(workspace->permutation); + + for(int i = 0; i < MAX_TEMP_MATRICES; i++) + gsl_matrix_free(workspace->tmp_mat[i]); + free(workspace->tmp_mat); + for(int i = 0; i < MAX_TEMP_VECTORS; i++) + gsl_vector_free(workspace->tmp_vec[i]); + free(workspace->tmp_vec); +} + +/************************************************** basic operations ********************************************************/ + +void invert(gsl_matrix *in, gsl_matrix *out, workspace_t *ws) +{ + int s; + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, in); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + gsl_linalg_LU_invert(tmp, ws->permutation, out); + + releaseTempMatrices(ws, 1); +} + +void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws) +{ + int s; + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, A); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + gsl_linalg_LU_solve(tmp, ws->permutation, b, result); + + releaseTempMatrices(ws, 1); +} + +void conjugate(gsl_matrix *in, gsl_matrix *conjugator, gsl_matrix *out, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + + invert(conjugator, out, ws); // use out to temporarily store inverse conjugator + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, in, out, 0.0, tmp); // in * conjugator^{-1} + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, conjugator, tmp, 0.0, out); + + releaseTempMatrices(ws, 1); +} + +void multiply(gsl_matrix *a, gsl_matrix *b, gsl_matrix *out) +{ + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, out); +} + +void multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, tmp); + gsl_matrix_memcpy(a, tmp); + releaseTempMatrices(ws, 1); +} + +void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, tmp); + gsl_matrix_memcpy(b, tmp); + releaseTempMatrices(ws, 1); +} + +void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...) +{ + va_list args; + va_start(args, n); + + gsl_matrix_set_identity(out); + + for(int i = 0; i < n; i++) { + gsl_matrix *cur = va_arg(args, gsl_matrix *); + multiply_right(out, cur, ws); + } + + va_end(args); +} + +void cartan_calc(gsl_matrix *g, double *mu, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, g); + gsl_linalg_SV_decomp(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)); + + releaseTempMatrices(ws, 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; + double result; + gsl_matrix *tmp = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, g); + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + result = gsl_linalg_LU_det(tmp, s); + + releaseTempMatrices(ws, 1); + return result; +} + + +int eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) +{ + gsl_matrix *g_ = getTempMatrix(ws); + int success = 0; + + gsl_matrix_memcpy(g_, g); + gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); + int r = gsl_eigen_nonsymmv(g_, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); + ERROR(r, "gsl_eigen_nonsymmv failed!\n"); + + gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); + + int real = 1; + for(int i = 0; i < ws->n; i++) + if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) != 0) + real = 0; + + if(!real) + goto eigenvectors_out; + + for(int i = 0; i < ws->n; i++) + for(int j = 0; j < ws->n; j++) + gsl_matrix_set(evec_real, i, j, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, i, j))); + + success = 1; + +eigenvectors_out: + releaseTempMatrices(ws, 1); + return success; +} + +// only fills in the real eigenvectors and returns their count +int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec_real, workspace_t *ws) +{ + gsl_matrix *g_ = getTempMatrix(ws); + + gsl_matrix_memcpy(g_, g); + gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); + int r = gsl_eigen_nonsymmv(g_, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); + ERROR(r, "gsl_eigen_nonsymmv failed!\n"); + + gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); + + int real = 0; + + for(int i = 0; i < ws->n; i++) { + if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, i)), 0) == 0) {// real + if(evec_real) { + for(int j = 0; j < ws->n; j++) + gsl_matrix_set(evec_real, j, real, GSL_REAL(gsl_matrix_complex_get(ws->evec_complex, j, i))); + } + real++; + } + } + + releaseTempMatrices(ws, 1); + return real; +} + +void eigensystem_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws) +{ + gsl_matrix *g_ = getTempMatrix(ws); + + gsl_matrix_memcpy(g_, g); + int r = gsl_eigen_symmv (g_, eval, evec, ws->work_symmv); + ERROR(r, "gsl_eigen_symmv failed!\n"); + + gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC); + + releaseTempMatrices(ws, 1); +} + +// returns number of positive directions and matrix transforming TO diagonal basis +int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws) +{ + gsl_matrix *A_ = getTempMatrix(ws); + + gsl_matrix_memcpy(A_, A); + int r = gsl_eigen_symmv (A_, ws->eval_real, cob, ws->work_symmv); + ERROR(r, "gsl_eigen_symmv failed!\n"); + + gsl_eigen_symmv_sort(ws->eval_real, cob, GSL_EIGEN_SORT_VAL_ASC); + + gsl_matrix_transpose(cob); + + int positive = 0; + + for(int i = 0; i < ws->n; i++) { + if(gsl_vector_get(ws->eval_real, i) > 0) + positive++; + + for(int j = 0; j < ws->n; j++) + *gsl_matrix_ptr(cob, i, j) *= sqrt(fabs(gsl_vector_get(ws->eval_real, i))); + } + + releaseTempMatrices(ws, 1); + return positive; +} + +// computes a matrix in SL(3, R) which projectively transforms (e1, e2, e3, e1+e2+e3) to the 4 given vectors +void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + gsl_vector *coeff = getTempVector(ws); + int s; + double det, scale; + + for(int i = 0; i < 3; i++) + for(int j = 0; j < 3; j++) + gsl_matrix_set(tmp, i, j, gsl_vector_get(vertices[j], i)); + + gsl_linalg_LU_decomp(tmp, ws->permutation, &s); + gsl_linalg_LU_solve(tmp, ws->permutation, vertices[3], coeff); + det = gsl_linalg_LU_det(tmp, s); + + for(int i = 0; i < 3; i++) + det *= gsl_vector_get(coeff, i); + scale = 1/cbrt(det); + + for(int i = 0; i < 3; i++) + for(int j = 0; j < 3; j++) + gsl_matrix_set(result, i, j, scale*gsl_vector_get(vertices[j], i)*gsl_vector_get(coeff, j)); + + releaseTempMatrices(ws, 1); + releaseTempVectors(ws, 1); +} + +void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws) +{ + gsl_matrix *tmp = getTempMatrix(ws); + gsl_matrix *rot_basis = getTempMatrix(ws); + + gsl_matrix_memcpy(tmp, rotation); + gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); + int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); + ERROR(r, "gsl_eigen_nonsymmv failed!\n"); + + double arg, minarg = 5; // greater than pi + int minidx; + for(int i = 0; i < 3; i++) { + arg = gsl_complex_arg(gsl_vector_complex_get(ws->eval_complex, i)); + if(abs(arg) < minarg) + { + minidx = i; + minarg = abs(arg); + } + } + ERROR(FCMP(minarg, 0.0) != 0, "rotation_frame() failed! No eigenvalue was 1.\n"); + + for(int i = 0; i < 3; i++) { + gsl_complex x = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+1)%3); + gsl_complex y = gsl_matrix_complex_get(ws->evec_complex, i, (minidx+2)%3); + gsl_complex z = gsl_matrix_complex_get(ws->evec_complex, i, minidx); + gsl_matrix_set(result, i, 0, GSL_REAL(x)+GSL_REAL(y)); + gsl_matrix_set(result, i, 1, GSL_IMAG(x)-GSL_IMAG(y)); + gsl_matrix_set(result, i, 2, GSL_REAL(z)); + } + + releaseTempMatrices(ws, 2); +} diff --git a/linalg.h b/linalg.h new file mode 100644 index 0000000..e077130 --- /dev/null +++ b/linalg.h @@ -0,0 +1,98 @@ +#ifndef LINALG_H +#define LINALG_H + +#include +#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) + +#define MAX_TEMP_MATRICES 10 +#define MAX_TEMP_VECTORS 10 + +typedef struct _workspace { + int n; + gsl_eigen_nonsymmv_workspace *work_nonsymmv; + gsl_eigen_symmv_workspace *work_symmv; + gsl_vector *work_sv; + gsl_vector_complex *eval_complex; + gsl_matrix_complex *evec_complex; + gsl_vector *eval_real; + gsl_matrix *evec_real; + gsl_permutation *permutation; + + gsl_matrix **tmp_mat; + int tmp_mat_used; + gsl_vector **tmp_vec; + int tmp_vec_used; +} workspace_t; + +workspace_t *workspace_alloc(int n); +void workspace_free(workspace_t *workspace); +void solve(gsl_matrix *A, gsl_vector *b, gsl_vector *result, workspace_t *ws); +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 multiply_right(gsl_matrix *a, gsl_matrix *b, workspace_t *ws); +void multiply_left(gsl_matrix *a, gsl_matrix *b, workspace_t *ws); +void multiply_many(workspace_t *ws, gsl_matrix *out, int n, ...); +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); +int jordan_calc(gsl_matrix *g, double *mu, workspace_t *ws); +double trace(gsl_matrix *g); +double determinant(gsl_matrix *g, workspace_t *ws); +int eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws); +int real_eigenvectors(gsl_matrix *g, gsl_matrix *evec, workspace_t *ws); +void eigenvectors_symm(gsl_matrix *g, gsl_vector *eval, gsl_matrix *evec, workspace_t *ws); +int diagonalize_symmetric_form(gsl_matrix *A, gsl_matrix *cob, workspace_t *ws); +void projective_frame(gsl_vector **vertices, gsl_matrix *result, workspace_t *ws); +void rotation_frame(gsl_matrix *rotation, gsl_matrix *result, workspace_t *ws); + +// matrix allocation stuff + +static gsl_matrix **getTempMatrices(workspace_t *ws, int n) +{ + ERROR(ws->tmp_mat_used + n > MAX_TEMP_MATRICES, "Ran out of temporary matrices. Consider increasing MAX_TEMP_MATRICES\n"); + int index = ws->tmp_mat_used; + ws->tmp_mat_used += n; + return ws->tmp_mat + index; +} + +static gsl_matrix *getTempMatrix(workspace_t *ws) +{ + return *getTempMatrices(ws, 1); +} + +static void releaseTempMatrices(workspace_t *ws, int n) +{ + ERROR(ws->tmp_mat_used - n < 0, "Released more matrices then in use\n"); + ws->tmp_mat_used -= n; +} + +static gsl_vector **getTempVectors(workspace_t *ws, int n) +{ + ERROR(ws->tmp_vec_used + n > MAX_TEMP_VECTORS, "Ran out of temporary vectors. Consider increasing MAX_TEMP_VECTORS\n"); + int index = ws->tmp_vec_used; + ws->tmp_vec_used += n; + return ws->tmp_vec + index; +} + +static gsl_vector *getTempVector(workspace_t *ws) +{ + return *getTempVectors(ws, 1); +} + +static void releaseTempVectors(workspace_t *ws, int n) +{ + ERROR(ws->tmp_vec_used - n < 0, "Released more vectors then in use\n"); + ws->tmp_vec_used -= n; +} + +#endif diff --git a/queue.h b/queue.h new file mode 100644 index 0000000..a37b2f2 --- /dev/null +++ b/queue.h @@ -0,0 +1,46 @@ +#ifndef QUEUE_H +#define QUEUE_H + +#include +#include + +#define QUEUE_SIZE 1500000 + +#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/singular_values.c b/singular_values.c new file mode 100644 index 0000000..44e9f76 --- /dev/null +++ b/singular_values.c @@ -0,0 +1,210 @@ +#include "triangle.h" +#include "linalg.h" + +#define MAX_ELEMENTS 14000 +//#define DRAW_PICTURE 1 + +void initialize_triangle_generators(workspace_t *ws, gsl_matrix **gen, double a1, double a2, double a3, double s, double t) +{ + gsl_matrix **tmp; + tmp = getTempMatrices(ws, 6); + + double rho1sqrt = sqrt(s*(s+2*cos(2*M_PI*a1))+1); + double rho2sqrt = sqrt(s*(s+2*cos(2*M_PI*a2))+1); + double rho3sqrt = sqrt(s*(s+2*cos(2*M_PI*a3))+1); + + for(int i = 0; i < 6; i++) + gsl_matrix_set_zero(tmp[i]); + + gsl_matrix_set(tmp[0], 0, 0, 1); + gsl_matrix_set(tmp[0], 1, 0, -rho3sqrt/t); + gsl_matrix_set(tmp[0], 2, 0, -rho2sqrt*t); + gsl_matrix_set(tmp[0], 1, 1, -1); + gsl_matrix_set(tmp[0], 2, 2, -1); + + gsl_matrix_set(tmp[1], 0, 0, -1); + gsl_matrix_set(tmp[1], 0, 1, -rho3sqrt*t); + gsl_matrix_set(tmp[1], 1, 1, 1); + gsl_matrix_set(tmp[1], 2, 1, -rho1sqrt/t); + gsl_matrix_set(tmp[1], 2, 2, -1); + + gsl_matrix_set(tmp[2], 0, 0, -1); + gsl_matrix_set(tmp[2], 1, 1, -1); + gsl_matrix_set(tmp[2], 0, 2, -rho2sqrt/t); + gsl_matrix_set(tmp[2], 1, 2, -rho1sqrt*t); + gsl_matrix_set(tmp[2], 2, 2, 1); + + gsl_matrix_set(tmp[3], 0, 0, 1); + gsl_matrix_set(tmp[3], 1, 1, s); + gsl_matrix_set(tmp[3], 2, 2, 1/s); + + gsl_matrix_set(tmp[4], 0, 0, 1/s); + gsl_matrix_set(tmp[4], 1, 1, 1); + gsl_matrix_set(tmp[4], 2, 2, s); + + gsl_matrix_set(tmp[5], 0, 0, s); + gsl_matrix_set(tmp[5], 1, 1, 1/s); + gsl_matrix_set(tmp[5], 2, 2, 1); + + multiply_many(ws, gen[0], 3, tmp[2], tmp[3], tmp[1]); + multiply_many(ws, gen[1], 3, tmp[0], tmp[4], tmp[2]); + multiply_many(ws, gen[2], 3, tmp[1], tmp[5], tmp[0]); + + invert(gen[0], gen[3], ws); + invert(gen[1], gen[4], ws); + invert(gen[2], gen[5], ws); + + releaseTempMatrices(ws, 6); +} + +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[6]; + gsl_matrix *tmp, *tmp2, *coxeter; + double mu[6], tr, trinv; + 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 < 6; 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); + +#ifdef DRAW_PICTURE + int width = 800; + int height = 800; + printf("P3\n%d %d\n255\n", width, height); + + for(int y = 0; y < height; y++) { + for(int x = 0; x < width; x++) { + double s = exp(((double)x/width-0.5)*4); + double t = exp(((double)y/height-0.5)*4); +#else + double s = atof(argv[1]); + double t = atof(argv[2]); +#endif + + // the real action + generate_triangle_group(group, MAX_ELEMENTS, 5, 6, 7); + initialize_triangle_generators(ws, gen, 1.0/5.0, 1.0/6.0, 1.0/7.0, s, t); + + gsl_matrix_set_identity(matrices[0]); + for(int i = 1; i < MAX_ELEMENTS; i++) { + if(group[i].length % 2 != 0) + continue; + + int parent = group[i].parent->id; + int grandparent = group[i].parent->parent->id; + int letter; + + if(group[parent].letter == 0 && group[i].letter == 1) + letter = 3; // p = ab + else if(group[parent].letter == 1 && group[i].letter == 2) + letter = 4; // q = bc + else if(group[parent].letter == 2 && group[i].letter == 0) + letter = 5; // r = ca + else if(group[parent].letter == 1 && group[i].letter == 0) + letter = 0; // p^{-1} = ba + else if(group[parent].letter == 2 && group[i].letter == 1) + letter = 1; // q^{-1} = cb + else if(group[parent].letter == 0 && group[i].letter == 2) + letter = 2; // r^{-1} = ac + + multiply(matrices[grandparent], gen[letter], matrices[i]); + } + + double excentricity; + double max_excentricity = 0; + int max_excentricity_index = 0; + + for(int i = 0; i < MAX_ELEMENTS; i++) { + if(group[i].length % 2 != 0) + continue; + + // jordan projection + gsl_matrix_memcpy(tmp, matrices[i]); + gsl_eigen_nonsymmv_params(0, ws->work_nonsymmv); + int r = gsl_eigen_nonsymmv(tmp, ws->eval_complex, ws->evec_complex, ws->work_nonsymmv); + ERROR(r, "gsl_eigen_nonsymmv failed!\n"); + gsl_eigen_nonsymmv_sort(ws->eval_complex, ws->evec_complex, GSL_EIGEN_SORT_ABS_DESC); + + int real_evs = 0; + + for(int j = 0; j < 3; j++) + if(FCMP(GSL_IMAG(gsl_vector_complex_get(ws->eval_complex, j)), 0) == 0) + real_evs++; + + if(real_evs != 3) + continue; + + mu[0] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 0)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1))); + mu[1] = fabs(GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 1)) / GSL_REAL(gsl_vector_complex_get(ws->eval_complex, 2))); + + excentricity = log(mu[1])/log(mu[0]); + if(excentricity > max_excentricity + 1e-3) { + max_excentricity = excentricity; + max_excentricity_index = i; + } + +#ifndef DRAW_PICTURE + printf("%d %f %f %f %f 0x0000ff %d %s\n", group[i].length, mu[0], mu[1], tr, trinv, i, print_word(&group[i], buf)); +#endif + } + +#ifdef DRAW_PICTURE + int color = (int)((max_excentricity-1)/(max_excentricity+4)*255); + if(color < 0) + color = 0; + if(color > 255) + color = 255; + printf("%d %d %d\n", color, color, color); + } + fprintf(stderr, "y = %d\n", y); + } +#else + fprintf(stderr, "max = %f %s\n", max_excentricity, print_word(&group[max_excentricity_index], buf)); +#endif + + // free stuff + for(int i = 0; i < MAX_ELEMENTS; i++) + gsl_matrix_free(matrices[i]); + for(int i = 0; i < 6; 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/singular_values.plt b/singular_values.plt new file mode 100644 index 0000000..d6775cc --- /dev/null +++ b/singular_values.plt @@ -0,0 +1,23 @@ +if(!exists("logt")) logt = log(1.5) +if(!exists("logs")) logs = log(1.0) + +file = sprintf("< ./singular_values %f %f", exp(logs), exp(logt)) +title = sprintf("s = %f, t = %f", exp(logs), exp(logt)) +# print title + +set zeroaxis +set size square +set xrange [0:25] +set yrange [0:25] +set grid +set parametric + +plot file using (log($2)):(log($3)) w p pt 7 ps 0.5 lc 1 t title +# plot file using 2:3 w p pt 7 ps 0.5 lc 1 t title + +pause mouse keypress +if(MOUSE_KEY == 60) logt=logt-0.02 +if(MOUSE_KEY == 62) logt=logt+0.02 +if(MOUSE_KEY == 44) logs=logs-0.02 +if(MOUSE_KEY == 46) logs=logs+0.02 +if(MOUSE_KEY != 113) reread diff --git a/triangle.c b/triangle.c new file mode 100644 index 0000000..770c7fe --- /dev/null +++ b/triangle.c @@ -0,0 +1,72 @@ +#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, *inv; + + int k[3] = {k1, k2, k3}; + + memset(group, 0, size*sizeof(groupelement_t)); + group[0].letter = -1; + 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; + } + } + + for(int i = 0; i < size; i++) { + cur = &group[i]; + inv = &group[0]; + while(cur->parent && inv) { + inv = inv->adj[cur->letter]; + cur = cur->parent; + } + group[i].inverse = inv; + } + + return n; +} diff --git a/triangle.h b/triangle.h new file mode 100644 index 0000000..076f6d4 --- /dev/null +++ b/triangle.h @@ -0,0 +1,21 @@ +#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; + struct _groupelement *inverse; + int letter; +} groupelement_t; + +int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3); + +#endif