triangle_reflection_complex/coxeter.c
2022-06-14 14:22:22 +02:00

204 lines
4.8 KiB
C

#include <math.h>
#include <malloc.h>
#include <string.h>
#include "coxeter.h"
#define LOOP(i,n) for(int i = 0; i < (n); i++)
#define MIN(x,y) ((x) > (y) ? (y) : (x))
#define MAX(x,y) ((x) > (y) ? (x) : (y))
group_t *coxeter_init_triangle(int p, int q, int r, int nmax)
{
int coxeter_matrix[9] = {
1, r, q,
r, 1, p,
q, p, 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);
}
int coxeter_snprint(char *str, int buflen, groupelement_t *g)
{
int n = MIN(g->length, buflen-1); // number of characters without null byte
str[n] = 0;
int i = n-1;
while(g->parent) {
str[i--] = 'a' + g->letter;
g = g->parent;
}
return n;
}