change Coxeter group algorithm
This commit is contained in:
Florian Stecker
parent
2568b4ef71
commit
9cb2e44867
16
Makefile
16
Makefile
|
|
@ -1,16 +1,16 @@
|
|||
HEADERS=triangle.h linalg.h queue.h mat.h
|
||||
HEADERS=linalg.h mat.h coxeter.h
|
||||
|
||||
#SPECIAL_OPTIONS=-O0 -g -D_DEBUG
|
||||
SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
|
||||
#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
|
||||
#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 mat.o
|
||||
gcc $(OPTIONS) -o singular_values triangle.o linalg.o singular_values.o mat.o -lm -lgsl -lcblas -lgmp -lmps -lpthread
|
||||
singular_values: singular_values.o coxeter.o linalg.o mat.o
|
||||
gcc $(OPTIONS) -o singular_values coxeter.o linalg.o singular_values.o mat.o -lm -lgsl -lcblas -lgmp -lmps -lpthread
|
||||
|
||||
singular_values.o: singular_values.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c singular_values.c
|
||||
|
|
@ -18,11 +18,11 @@ singular_values.o: singular_values.c $(HEADERS)
|
|||
linalg.o: linalg.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c linalg.c
|
||||
|
||||
triangle.o: triangle.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c triangle.c
|
||||
coxeter.o: coxeter.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c coxeter.c
|
||||
|
||||
mat.o: mat.c $(HEADERS)
|
||||
gcc $(OPTIONS) -c mat.c
|
||||
|
||||
clean:
|
||||
rm -f singular_values triangle.o linalg.o singular_values.o mat.o
|
||||
rm -f singular_values coxeter.o linalg.o singular_values.o mat.o
|
||||
|
|
|
|||
|
|
@ -0,0 +1,187 @@
|
|||
#include <math.h>
|
||||
#include <malloc.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "coxeter.h"
|
||||
|
||||
#define LOOP(i,n) for(int i = 0; i < (n); i++)
|
||||
|
||||
group_t *coxeter_init_triangle(int p, int q, int r, int nmax)
|
||||
{
|
||||
int coxeter_matrix[9] = {
|
||||
1, p, r,
|
||||
p, 1, q,
|
||||
r, q, 1};
|
||||
|
||||
return coxeter_init(3, coxeter_matrix, nmax);
|
||||
}
|
||||
|
||||
group_t *coxeter_init(int rank, int *coxeter_matrix, int nmax)
|
||||
{
|
||||
int n;
|
||||
group_t *group;
|
||||
double *schlaefli_matrix;
|
||||
double *vectors;
|
||||
double *cur_vec;
|
||||
groupelement_t *cur, *cur_inv;
|
||||
int word[100]; // reasonable estimate?
|
||||
|
||||
// allocate
|
||||
group = malloc(sizeof(group_t));
|
||||
group->coxeter_matrix = malloc(rank*rank*sizeof(int));
|
||||
schlaefli_matrix = malloc(rank*rank*sizeof(double));
|
||||
vectors = malloc(rank*nmax*sizeof(double));
|
||||
cur_vec = malloc(rank*sizeof(double));
|
||||
group->elements = malloc(nmax*sizeof(groupelement_t));
|
||||
group->lists = malloc(2*nmax*rank*sizeof(groupelement_t*));
|
||||
memset(group->lists, 0, 2*nmax*rank*sizeof(groupelement_t*));
|
||||
LOOP(i, nmax) group->elements[i].left = group->lists + i*rank;
|
||||
LOOP(i, nmax) group->elements[i].right = group->lists + (nmax+i)*rank;
|
||||
LOOP(i, nmax) group->elements[i].letter = -1;
|
||||
LOOP(i, nmax) group->elements[i].id = i;
|
||||
|
||||
// copy coxeter matrix
|
||||
group->rank = rank;
|
||||
memcpy(group->coxeter_matrix, coxeter_matrix, rank*rank*sizeof(int));
|
||||
|
||||
// generate Schläfli matrix
|
||||
LOOP(i,rank) LOOP(j,rank) {
|
||||
if(group->coxeter_matrix[i*rank+j] == -1)
|
||||
schlaefli_matrix[i*rank+j] = -2;
|
||||
else
|
||||
schlaefli_matrix[i*rank+j] = -2*cos(M_PI/group->coxeter_matrix[i*rank+j]);
|
||||
}
|
||||
|
||||
// identity element
|
||||
group->elements[0].length = 0;
|
||||
group->elements[0].letter = -1;
|
||||
LOOP(i, rank) vectors[i] = 1.0;
|
||||
n = 1;
|
||||
|
||||
// elements and left multiplication
|
||||
for(int i = 0; n < nmax && i < n; i++) {
|
||||
LOOP(j, rank) {
|
||||
if(n >= nmax)
|
||||
break;
|
||||
if(vectors[i*rank+j] < 0) // this generator decreases length
|
||||
continue;
|
||||
if(group->elements[i].left[j] != 0) // we already added this element
|
||||
continue;
|
||||
|
||||
LOOP(k, rank) {
|
||||
if(k == j)
|
||||
vectors[n*rank+k] = -vectors[i*rank+k];
|
||||
else
|
||||
vectors[n*rank+k] = vectors[i*rank+k] - vectors[i*rank+j]*schlaefli_matrix[k*rank+j];
|
||||
}
|
||||
|
||||
group->elements[n].length = group->elements[i].length + 1;
|
||||
|
||||
// if s_k * w is shorter than w (the new element), find out what it is and update "left"
|
||||
LOOP(k, rank) {
|
||||
if(vectors[n*rank+k] > 0)
|
||||
continue;
|
||||
|
||||
if(group->elements[n].letter == -1)
|
||||
group->elements[n].letter = k;
|
||||
|
||||
// get w
|
||||
LOOP(l, rank) cur_vec[l] = vectors[n*rank+l];
|
||||
|
||||
// apply s_k
|
||||
LOOP(l, rank)
|
||||
if(l != k)
|
||||
cur_vec[l] -= cur_vec[k]*schlaefli_matrix[l*rank+k];
|
||||
cur_vec[k] = -cur_vec[k];
|
||||
|
||||
// find a reduced word for s_k w
|
||||
LOOP(m, group->elements[i].length) { // s_k w should have same length as element i
|
||||
int p;
|
||||
|
||||
// find a generator s_p decreasing the word length
|
||||
for(p = 0; p < rank; p++)
|
||||
if(cur_vec[p] < 0)
|
||||
break;
|
||||
|
||||
if(p == rank) // this should not happen
|
||||
fprintf(stderr, "Uh oh!\n");
|
||||
|
||||
word[m] = p;
|
||||
|
||||
// apply s_p
|
||||
LOOP(l, rank)
|
||||
if(l != p)
|
||||
cur_vec[l] -= cur_vec[p]*schlaefli_matrix[l*rank+p];
|
||||
cur_vec[p] = -cur_vec[p];
|
||||
}
|
||||
|
||||
// find the element corresponding to the word
|
||||
groupelement_t *cur = &group->elements[0];
|
||||
LOOP(m, group->elements[i].length) {
|
||||
cur = cur->left[word[group->elements[i].length - 1 - m]];
|
||||
}
|
||||
|
||||
cur->left[k] = &group->elements[n];
|
||||
group->elements[n].left[k] = cur;
|
||||
}
|
||||
|
||||
n++;
|
||||
}
|
||||
}
|
||||
|
||||
group->size = n;
|
||||
|
||||
// parent
|
||||
LOOP(i, n) {
|
||||
if(i == 0)
|
||||
group->elements[i].parent = NULL;
|
||||
else
|
||||
group->elements[i].parent = group->elements[i].left[group->elements[i].letter];
|
||||
}
|
||||
|
||||
// right multiplication
|
||||
LOOP(i, n) {
|
||||
LOOP(j, rank) {
|
||||
if(group->elements[i].right[j])
|
||||
continue;
|
||||
|
||||
int k = 0;
|
||||
for(groupelement_t *cur = &group->elements[i]; cur->parent; cur = cur->parent)
|
||||
word[k++] = cur->letter;
|
||||
cur = group->elements[0].left[j];
|
||||
for(int k = group->elements[i].length - 1; k >= 0; k--)
|
||||
if(cur)
|
||||
cur = cur->left[word[k]];
|
||||
|
||||
if(cur) {
|
||||
group->elements[i].right[j] = cur;
|
||||
cur->right[j] = &group->elements[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// inverse
|
||||
LOOP(i, n) {
|
||||
cur_inv = &group->elements[0];
|
||||
for(groupelement_t *cur = &group->elements[i]; cur->parent; cur = cur->parent) {
|
||||
if(cur_inv == NULL)
|
||||
break;
|
||||
cur_inv = cur_inv->left[cur->letter];
|
||||
}
|
||||
group->elements[i].inverse = cur_inv;
|
||||
}
|
||||
|
||||
// free
|
||||
free(schlaefli_matrix);
|
||||
free(vectors);
|
||||
|
||||
return group;
|
||||
}
|
||||
|
||||
void coxeter_clear(group_t *g)
|
||||
{
|
||||
free(g->coxeter_matrix);
|
||||
free(g->elements);
|
||||
free(g->lists);
|
||||
free(g);
|
||||
}
|
||||
|
|
@ -0,0 +1,29 @@
|
|||
#ifndef COXETER_H
|
||||
#define COXETER_H
|
||||
|
||||
#include <inttypes.h>
|
||||
|
||||
typedef struct _groupelement {
|
||||
int id;
|
||||
int letter;
|
||||
struct _groupelement *parent;
|
||||
|
||||
int length;
|
||||
struct _groupelement *inverse;
|
||||
struct _groupelement **left;
|
||||
struct _groupelement **right;
|
||||
} groupelement_t;
|
||||
|
||||
typedef struct _group {
|
||||
int rank;
|
||||
int *coxeter_matrix;
|
||||
int size;
|
||||
groupelement_t *elements;
|
||||
groupelement_t **lists;
|
||||
} group_t;
|
||||
|
||||
group_t *coxeter_init(int rank, int *coxeter_matrix, int nmax);
|
||||
group_t *coxeter_init_triangle(int p, int q, int r, int nmax);
|
||||
void coxeter_clear(group_t *g);
|
||||
|
||||
#endif
|
||||
46
queue.h
46
queue.h
|
|
@ -1,46 +0,0 @@
|
|||
#ifndef QUEUE_H
|
||||
#define QUEUE_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#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
|
||||
|
|
@ -1,16 +1,32 @@
|
|||
#include "triangle.h"
|
||||
#include "coxeter.h"
|
||||
#include "linalg.h"
|
||||
#include "mat.h"
|
||||
|
||||
#include <gsl/gsl_poly.h>
|
||||
#include <mps/mps.h>
|
||||
|
||||
//#define MAX_ELEMENTS 2800000
|
||||
//#define MAX_ELEMENTS 720000
|
||||
#define MAX_ELEMENTS 1000
|
||||
//#define DRAW_PICTURE 1
|
||||
|
||||
#define SWAP(t,x,y) do { t _tmp = (x); (x) = (y); (y) = _tmp; } while (0);
|
||||
//#define DEBUG(msg, ...) fprintf(stderr, msg, ##__VA_ARGS__)
|
||||
#define DEBUG(msg, ...)
|
||||
|
||||
struct result {
|
||||
mpq_t tr;
|
||||
mpq_t trinv;
|
||||
};
|
||||
|
||||
static int compare_result(const void *a_, const void *b_)
|
||||
{
|
||||
int d = 0;
|
||||
|
||||
struct result **a = (struct result **)a_;
|
||||
struct result **b = (struct result **)b_;
|
||||
|
||||
d = mpq_cmp((*a)->tr,(*b)->tr);
|
||||
if(d == 0)
|
||||
mpq_cmp((*a)->trinv, (*b)->trinv);
|
||||
|
||||
return d;
|
||||
}
|
||||
|
||||
int solve_characteristic_polynomial(mps_context *solv, mpq_t tr, mpq_t trinv, double *eigenvalues)
|
||||
{
|
||||
|
|
@ -213,7 +229,7 @@ char *print_word(groupelement_t *g, char *str)
|
|||
return str;
|
||||
}
|
||||
|
||||
void enumerate(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t)
|
||||
void enumerate(group_t *group, mat *matrices, mpq_t s, mpq_t t)
|
||||
{
|
||||
mat_workspace *ws;
|
||||
mat tmp;
|
||||
|
|
@ -229,27 +245,27 @@ void enumerate(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t)
|
|||
initialize_triangle_generators(ws, gen, s, t);
|
||||
|
||||
mat_identity(matrices[0]);
|
||||
for(int i = 1; i < MAX_ELEMENTS; i++) {
|
||||
if(group[i].length % 2 != 0)
|
||||
for(int i = 1; i < group->size; i++) {
|
||||
if(group->elements[i].length % 2 != 0)
|
||||
continue;
|
||||
if(!group[i].inverse)
|
||||
if(!group->elements[i].inverse)
|
||||
continue;
|
||||
|
||||
int parent = group[i].parent->id;
|
||||
int grandparent = group[i].parent->parent->id;
|
||||
int parent = group->elements[i].parent->id;
|
||||
int grandparent = group->elements[i].parent->parent->id;
|
||||
int letter;
|
||||
|
||||
if(group[parent].letter == 1 && group[i].letter == 2)
|
||||
if(group->elements[parent].letter == 1 && group->elements[i].letter == 2)
|
||||
letter = 0; // p = bc
|
||||
else if(group[parent].letter == 2 && group[i].letter == 0)
|
||||
else if(group->elements[parent].letter == 2 && group->elements[i].letter == 0)
|
||||
letter = 1; // q = ca
|
||||
else if(group[parent].letter == 0 && group[i].letter == 1)
|
||||
else if(group->elements[parent].letter == 0 && group->elements[i].letter == 1)
|
||||
letter = 2; // r = ab
|
||||
if(group[parent].letter == 2 && group[i].letter == 1)
|
||||
if(group->elements[parent].letter == 2 && group->elements[i].letter == 1)
|
||||
letter = 3; // p^{-1} = cb
|
||||
else if(group[parent].letter == 0 && group[i].letter == 2)
|
||||
else if(group->elements[parent].letter == 0 && group->elements[i].letter == 2)
|
||||
letter = 4; // q^{-1} = ac
|
||||
else if(group[parent].letter == 1 && group[i].letter == 0)
|
||||
else if(group->elements[parent].letter == 1 && group->elements[i].letter == 0)
|
||||
letter = 5; // r^{-1} = ba
|
||||
|
||||
mat_multiply(ws, matrices[i], matrices[grandparent], gen[letter]);
|
||||
|
|
@ -262,7 +278,7 @@ void enumerate(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t)
|
|||
mat_workspace_clear(ws);
|
||||
}
|
||||
|
||||
void output_invariants(groupelement_t *group, mat *matrices, mpq_t s, mpq_t q, mps_context *solver)
|
||||
void output_invariants(group_t *group, mat *matrices, mpq_t s, mpq_t q, mps_context *solver)
|
||||
{
|
||||
mpq_t tr, trinv;
|
||||
char buf[100];
|
||||
|
|
@ -271,12 +287,12 @@ void output_invariants(groupelement_t *group, mat *matrices, mpq_t s, mpq_t q, m
|
|||
|
||||
mpq_inits(tr, trinv, NULL);
|
||||
|
||||
for(int i = 0; i < MAX_ELEMENTS; i++) {
|
||||
if(group[i].length % 2 != 0 || !group[i].inverse)
|
||||
for(int i = 0; i < group->size; i++) {
|
||||
if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse)
|
||||
continue;
|
||||
|
||||
mat_trace(tr, matrices[i]);
|
||||
mat_trace(trinv, matrices[group[i].inverse->id]);
|
||||
mat_trace(trinv, matrices[group->elements[i].inverse->id]);
|
||||
|
||||
retval = solve_characteristic_polynomial(solver, tr, trinv, evs);
|
||||
if(retval == 1) {
|
||||
|
|
@ -293,12 +309,13 @@ void output_invariants(groupelement_t *group, mat *matrices, mpq_t s, mpq_t q, m
|
|||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
|
||||
gmp_printf("%d %d %s %Qd %Qd %f %f\n", i, group[i].length, print_word(&group[i], buf), tr, trinv, log(evs[0]), -log(evs[2]));
|
||||
gmp_printf("%d %d %s %Qd %Qd %f %f\n", i, group->elements[i].length, print_word(&group->elements[i], buf), tr, trinv, log(evs[0]), -log(evs[2]));
|
||||
}
|
||||
|
||||
mpq_clears(tr, trinv, NULL);
|
||||
}
|
||||
|
||||
/*
|
||||
double max_slope(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t, int *index)
|
||||
{
|
||||
double max = 0;
|
||||
|
|
@ -328,37 +345,58 @@ double max_slope(groupelement_t *group, mat *matrices, mpq_t s, mpq_t t, int *in
|
|||
|
||||
return max;
|
||||
}
|
||||
*/
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
mpq_t s, q, t, tmp;
|
||||
double sapprox, tapprox, qapprox, tqfactor;
|
||||
mat *matrices;
|
||||
groupelement_t *group;
|
||||
group_t *group;
|
||||
int index;
|
||||
mps_context *solver;
|
||||
int acc = 100;
|
||||
int n;
|
||||
int retval;
|
||||
double evs[3];
|
||||
double max_slope;
|
||||
int nmax;
|
||||
|
||||
struct result *invariants;
|
||||
struct result **distinct_invariants;
|
||||
|
||||
nmax = atoi(argv[1]);
|
||||
|
||||
DEBUG("Allocate\n");
|
||||
|
||||
mpq_inits(s, q, t, tmp, NULL);
|
||||
group = malloc(MAX_ELEMENTS*sizeof(groupelement_t));
|
||||
matrices = malloc(MAX_ELEMENTS*sizeof(mat));
|
||||
for(int i = 0; i < MAX_ELEMENTS; i++)
|
||||
matrices = malloc(nmax*sizeof(mat));
|
||||
for(int i = 0; i < nmax; i++)
|
||||
mat_init(matrices[i], 3);
|
||||
invariants = malloc(nmax*sizeof(struct result));
|
||||
distinct_invariants = malloc(nmax*sizeof(struct result));
|
||||
for(int i = 0; i < nmax; i++) {
|
||||
mpq_init(invariants[i].tr);
|
||||
mpq_init(invariants[i].trinv);
|
||||
distinct_invariants[i] = &invariants[i];
|
||||
}
|
||||
|
||||
solver = mps_context_new();
|
||||
mps_context_set_output_prec(solver, 20); // relative precision
|
||||
mps_context_set_output_goal(solver, MPS_OUTPUT_GOAL_APPROXIMATE);
|
||||
|
||||
int acc = 100;
|
||||
DEBUG("Approximate parameters\n");
|
||||
|
||||
sapprox = atof(argv[1]);
|
||||
tapprox = atof(argv[2]);
|
||||
// get approximate s and q values
|
||||
sapprox = atof(argv[2]);
|
||||
tapprox = atof(argv[3]);
|
||||
tqfactor = pow((sapprox*sapprox-sapprox+1)*(sapprox*sapprox-sapprox+1)*(sapprox*sapprox+1), 1/6.0);
|
||||
qapprox = tapprox/tqfactor;
|
||||
|
||||
for(int i = 0; ; i++) {
|
||||
continued_fraction_approximation(tmp, sapprox, i);
|
||||
// gmp_fprintf(stdout, "%Qd\n", tmp);
|
||||
if(fabs(mpq_get_d(t)-sapprox) < 1e-10 || (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0))
|
||||
if(fabs(mpq_get_d(t)-sapprox) < 1e-10
|
||||
|| (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0))
|
||||
break;
|
||||
mpq_set(s, tmp);
|
||||
}
|
||||
|
|
@ -366,8 +404,8 @@ int main(int argc, char *argv[])
|
|||
|
||||
for(int i = 0; ; i++) {
|
||||
continued_fraction_approximation(tmp, qapprox, i);
|
||||
// gmp_fprintf(stdout, "%Qd\n", tmp);
|
||||
if(fabs(mpq_get_d(t)-qapprox) < 1e-10 || (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0))
|
||||
if(fabs(mpq_get_d(t)-qapprox) < 1e-10
|
||||
|| (mpz_cmpabs_ui(mpq_numref(tmp),acc) > 0 && mpz_cmpabs_ui(mpq_denref(tmp),acc) > 0))
|
||||
break;
|
||||
mpq_set(q, tmp);
|
||||
}
|
||||
|
|
@ -375,23 +413,83 @@ int main(int argc, char *argv[])
|
|||
|
||||
tqfactor = pow((mpq_get_d(s)*mpq_get_d(s)-mpq_get_d(s)+1)*(mpq_get_d(s)*mpq_get_d(s)-mpq_get_d(s)+1)*(mpq_get_d(s)*mpq_get_d(s)+1), 1/6.0);
|
||||
|
||||
gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor);
|
||||
// gmp_fprintf(stdout, "\"s = %Qd = %.3f, q = %Qd, t = %.3f\"\n", s, mpq_get_d(s), q, mpq_get_d(q)*tqfactor);
|
||||
|
||||
generate_triangle_group(group, MAX_ELEMENTS, 3, 3, 4);
|
||||
// group
|
||||
DEBUG("Generate group\n");
|
||||
group = coxeter_init_triangle(3, 3, 4, nmax);
|
||||
|
||||
DEBUG("Compute matrices\n");
|
||||
enumerate(group, matrices, s, q);
|
||||
|
||||
DEBUG("Compute traces\n");
|
||||
for(int i = 0; i < nmax; i++) {
|
||||
if(group->elements[i].length % 2 != 0 || !group->elements[i].inverse)
|
||||
continue;
|
||||
|
||||
mat_trace(invariants[i].tr, matrices[i]);
|
||||
mat_trace(invariants[i].trinv, matrices[group->elements[i].inverse->id]);
|
||||
}
|
||||
|
||||
DEBUG("Get unique traces\n");
|
||||
qsort(distinct_invariants, nmax, sizeof(struct result*), compare_result);
|
||||
|
||||
n = 0;
|
||||
for(int i = 0; i < nmax; i++) {
|
||||
if(i == 0 || compare_result(&distinct_invariants[i], &distinct_invariants[n-1]) != 0)
|
||||
distinct_invariants[n++] = distinct_invariants[i];
|
||||
}
|
||||
|
||||
max_slope = 0;
|
||||
|
||||
DEBUG("Solve characteristic polynomials\n");
|
||||
for(int i = 0; i < n; i++) {
|
||||
retval = solve_characteristic_polynomial(solver, distinct_invariants[i]->tr, distinct_invariants[i]->trinv, evs);
|
||||
if(retval == 1) {
|
||||
fprintf(stderr, "Error! Could not solve polynomial.\n");
|
||||
continue;
|
||||
} else if(retval == 2) {
|
||||
continue;
|
||||
}
|
||||
|
||||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
if(fabs(evs[1]) < fabs(evs[2]))
|
||||
SWAP(double, evs[1], evs[2]);
|
||||
if(fabs(evs[0]) < fabs(evs[1]))
|
||||
SWAP(double, evs[0], evs[1]);
|
||||
|
||||
double x = log(evs[0]);
|
||||
double y = -log(evs[2]);
|
||||
|
||||
if(y/x > max_slope && (x > 0.1 || y > 0.1))
|
||||
max_slope = y/x;
|
||||
|
||||
gmp_printf("%Qd %Qd %f %f %f\n", distinct_invariants[i]->tr, distinct_invariants[i]->trinv, x, y, y/x);
|
||||
}
|
||||
|
||||
// printf("%.3f %.3f %f\n", mpq_get_d(s), mpq_get_d(q)*tqfactor, max_slope);
|
||||
|
||||
// output_invariants(group, matrices, s, q, solver);
|
||||
|
||||
// for(int i = 0; i < 10; i++) {
|
||||
// mpq_set_ui(t,100+i,100);
|
||||
// mpq_canonicalize(t);
|
||||
|
||||
enumerate(group, matrices, s, q);
|
||||
//printf("%f %f\n", mpq_get_d(t), max_slope(group, matrices, s, t, &index));
|
||||
output_invariants(group, matrices, s, q, solver);
|
||||
// }
|
||||
|
||||
for(int i = 0; i < MAX_ELEMENTS; i++)
|
||||
DEBUG("Clean up\n");
|
||||
for(int i = 0; i < nmax; i++) {
|
||||
mpq_clear(invariants[i].tr);
|
||||
mpq_clear(invariants[i].trinv);
|
||||
}
|
||||
free(invariants);
|
||||
free(distinct_invariants);
|
||||
for(int i = 0; i < nmax; i++)
|
||||
mat_clear(matrices[i]);
|
||||
free(matrices);
|
||||
free(group);
|
||||
coxeter_clear(group);
|
||||
mpq_clears(s, q, t, tmp, NULL);
|
||||
mps_context_free(solver);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
if(!exists("logt")) logt = log(1.78)
|
||||
if(!exists("logs")) logs = log(0.5)
|
||||
|
||||
file = sprintf("< ./singular_values %f %f", exp(logs), exp(logt))
|
||||
file = sprintf("< ./singular_values 1000 %f %f", exp(logs), exp(logt))
|
||||
|
||||
set zeroaxis
|
||||
set samples 1000
|
||||
|
|
@ -24,7 +24,7 @@ trinv(a,b) = exp(-a) + exp(a-b) + exp(b)
|
|||
# log(tr(t,t*2)),log(trinv(t,2*t)) w l lw 2 t "", \
|
||||
# log(tr(t,t/2)),log(trinv(t,t/2)) w l lw 2 t ""
|
||||
|
||||
plot file using 6:7 w p pt 7 ps 0.5 lc 1 t columnheader, \
|
||||
plot file using 3:4 w p pt 7 ps 0.5 lc 1 t columnheader, \
|
||||
t,2*t w l lw 2 t "", \
|
||||
t,t/2 w l lw 2 t ""
|
||||
|
||||
|
|
|
|||
72
triangle.c
72
triangle.c
|
|
@ -1,72 +0,0 @@
|
|||
#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, *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;
|
||||
}
|
||||
21
triangle.h
21
triangle.h
|
|
@ -1,21 +0,0 @@
|
|||
#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;
|
||||
struct _groupelement *inverse;
|
||||
int letter;
|
||||
} groupelement_t;
|
||||
|
||||
int generate_triangle_group(groupelement_t *group, int size, int k1, int k2, int k3);
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue