Thickenings
This commit is contained in:
Florian Stecker
commit
53811ef32a
14
Makefile
Normal file
14
Makefile
Normal file
@ -0,0 +1,14 @@
|
||||
|
||||
all: thickenings
|
||||
|
||||
thickenings: thickenings.o coxeter.o
|
||||
gcc -g -o thickenings thickenings.o coxeter.o -lgsl -lcblas
|
||||
|
||||
thickenings.o: thickenings.c coxeter.h
|
||||
gcc -g -c thickenings.c -std=gnu99
|
||||
|
||||
coxeter.o: coxeter.c coxeter.h
|
||||
gcc -g -c coxeter.c -std=gnu99
|
||||
|
||||
clean:
|
||||
rm -f thickenings thickenings.o coxeter.o
|
||||
242
coxeter.c
Normal file
242
coxeter.c
Normal file
@ -0,0 +1,242 @@
|
||||
#include "coxeter.h"
|
||||
#include "queue.h"
|
||||
|
||||
#include <memory.h>
|
||||
#include <gsl/gsl_linalg.h>
|
||||
#include <gsl/gsl_blas.h>
|
||||
|
||||
static void applyGenerator(int i, double *x, double *y, int rank, double **schlaefliMatrix)
|
||||
{
|
||||
memcpy(y, x, rank*sizeof(double));
|
||||
for(int j = 0; j < rank; j++)
|
||||
y[i] -= schlaefliMatrix[i][j]*x[j];
|
||||
}
|
||||
|
||||
static int check_equal_sector(gsl_matrix *x1, gsl_matrix *x2, gsl_matrix *yinv, gsl_matrix *schlaefli, gsl_vector *starting, gsl_vector *i1, gsl_vector *i2)
|
||||
{
|
||||
double scalarProduct;
|
||||
int rank = x1->size1;
|
||||
|
||||
// calculate <y^{-1} * x1 * x2 * s, e_i> for all i, or identically the components of schlaefli*y^{-1}*x*s
|
||||
gsl_blas_dgemv(CblasNoTrans, 1.0, x2, starting, 0.0, i1);
|
||||
gsl_blas_dgemv(CblasNoTrans, 1.0, x1, i1, 0.0, i2);
|
||||
gsl_blas_dgemv(CblasNoTrans, 1.0, yinv, i2, 0.0, i1);
|
||||
gsl_blas_dgemv(CblasNoTrans, 1.0, schlaefli, i1, 0.0, i2);
|
||||
for(int i = 0; i < rank; i++)
|
||||
if(gsl_vector_get(i2, i) < 0)
|
||||
return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
static void invert(gsl_matrix *x, gsl_matrix *xinv)
|
||||
{
|
||||
int rank = x->size1;
|
||||
gsl_matrix *xLU = gsl_matrix_alloc(rank, rank);
|
||||
gsl_permutation *p = gsl_permutation_alloc(rank);
|
||||
int s;
|
||||
gsl_matrix_memcpy(xLU, x);
|
||||
gsl_linalg_LU_decomp(xLU, p, &s);
|
||||
gsl_linalg_LU_invert(xLU, p, xinv);
|
||||
gsl_permutation_free(p);
|
||||
gsl_matrix_free(xLU);
|
||||
}
|
||||
|
||||
static void generate_schlaefli_matrix(semisimple_type_t type, gsl_matrix *mat)
|
||||
{
|
||||
gsl_matrix_set_zero(mat);
|
||||
int offset = 0;
|
||||
|
||||
for(int k = 0; k < type.n; k++) {
|
||||
if(type.factors[k].series == 'A') {
|
||||
for(int i = 0; i < type.factors[k].rank; i++)
|
||||
for(int j = 0; j < type.factors[k].rank; j++)
|
||||
if(i == j)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, 2.0);
|
||||
else if(i - j == 1 || i - j == -1)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, -1.0);
|
||||
} else if(type.factors[k].series == 'B' || type.factors[k].series == 'C') {
|
||||
for(int i = 0; i < type.factors[k].rank; i++)
|
||||
for(int j = 0; j < type.factors[k].rank; j++)
|
||||
if(i == j)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, 2.0);
|
||||
else if(i == 0 && j == 1 || i == 1 && j == 0)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, -sqrt(2.0));
|
||||
else if(i - j == 1 || i - j == -1)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, -1.0);
|
||||
} else if(type.factors[k].series == 'D') {
|
||||
for(int i = 0; i < type.factors[k].rank; i++)
|
||||
for(int j = 0; j < type.factors[k].rank; j++)
|
||||
if(i == j)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, 2.0);
|
||||
else if(i > 0 && j > 0 && (i - j == 1 || i - j == -1) || i == 0 && j == 2 || i == 2 && j == 0)
|
||||
gsl_matrix_set(mat, offset + i, offset + j, -1.0);
|
||||
} else if(type.factors[k].series == 'G') {
|
||||
ERROR(type.factors[k].rank != 2, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
|
||||
gsl_matrix_set(mat, offset + 0, offset + 0, 2.0);
|
||||
gsl_matrix_set(mat, offset + 0, offset + 1, -sqrt(3.0));
|
||||
gsl_matrix_set(mat, offset + 1, offset + 0, -sqrt(3.0));
|
||||
gsl_matrix_set(mat, offset + 1, offset + 1, 2.0);
|
||||
} else {
|
||||
ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
|
||||
}
|
||||
|
||||
offset += type.factors[k].rank;
|
||||
}
|
||||
}
|
||||
|
||||
int coxeter_rank(semisimple_type_t type)
|
||||
{
|
||||
int rank = 0;
|
||||
for(int i = 0; i < type.n; i++)
|
||||
rank += type.factors[i].rank;
|
||||
return rank;
|
||||
}
|
||||
|
||||
int coxeter_order(semisimple_type_t type)
|
||||
{
|
||||
int order = 1;
|
||||
for(int i = 0; i < type.n; i++) {
|
||||
if(type.factors[i].series == 'A') { // (rank+1)!
|
||||
for(int j = 1; j <= type.factors[i].rank + 1; j++)
|
||||
order *= j;
|
||||
} else if(type.factors[i].series == 'B' || type.factors[i].series == 'C') { // 2^rank * rank!
|
||||
for(int j = 1; j <= type.factors[i].rank; j++)
|
||||
order *= 2*j;
|
||||
} else if(type.factors[i].series == 'D') { // 2^(rank-1) * rank!
|
||||
for(int j = 2; j <= type.factors[i].rank; j++)
|
||||
order *= 2*j;
|
||||
} else if(type.factors[i].series == 'E') {
|
||||
if(type.factors[i].rank == 6)
|
||||
order *= 51840;
|
||||
else if(type.factors[i].rank == 7)
|
||||
order *= 2903040;
|
||||
else if(type.factors[i].rank == 8)
|
||||
order *= 696729600;
|
||||
else
|
||||
ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
|
||||
} else if(type.factors[i].series == 'F') {
|
||||
ERROR(type.factors[i].rank != 4, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
|
||||
order *= 1152;
|
||||
} else if(type.factors[i].series == 'G') {
|
||||
ERROR(type.factors[i].rank != 2, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
|
||||
order *= 12;
|
||||
} else if(type.factors[i].series == 'H') {
|
||||
if(type.factors[i].rank == 2)
|
||||
order *= 10;
|
||||
else if(type.factors[i].rank == 3)
|
||||
order *= 120;
|
||||
else if(type.factors[i].rank == 4)
|
||||
order *= 14400;
|
||||
else
|
||||
ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
|
||||
} else {
|
||||
ERROR(1, "A Coxeter group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
|
||||
}
|
||||
}
|
||||
return order;
|
||||
}
|
||||
|
||||
static void generate_starting_vector(int rank, gsl_matrix *schlaefli, gsl_vector *result)
|
||||
{
|
||||
gsl_matrix *schlaefliLU = gsl_matrix_alloc(rank, rank);
|
||||
gsl_vector *diagonal = gsl_vector_alloc(rank);
|
||||
gsl_permutation *p = gsl_permutation_alloc(rank);
|
||||
int s;
|
||||
for(int i = 0; i < rank; i++)
|
||||
gsl_vector_set(diagonal, i, 1.0);
|
||||
gsl_matrix_memcpy(schlaefliLU, schlaefli);
|
||||
gsl_linalg_LU_decomp(schlaefliLU, p, &s);
|
||||
gsl_linalg_LU_solve(schlaefliLU, p, diagonal, result);
|
||||
gsl_permutation_free(p);
|
||||
gsl_vector_free(diagonal);
|
||||
gsl_matrix_free(schlaefliLU);
|
||||
}
|
||||
|
||||
void generate_coxeter_graph(semisimple_type_t type, int *result)
|
||||
{
|
||||
int rank = coxeter_rank(type);
|
||||
int order = coxeter_order(type);
|
||||
int element_count;
|
||||
gsl_matrix *schlaefliMatrix = gsl_matrix_alloc(rank, rank);
|
||||
gsl_vector *startingVector = gsl_vector_alloc(rank);
|
||||
gsl_matrix *generators = (gsl_matrix*)malloc(rank*sizeof(gsl_matrix));
|
||||
double *generators_data = (double*)malloc(rank*rank*rank*sizeof(double));
|
||||
gsl_matrix *elements = (gsl_matrix*)malloc(order*sizeof(gsl_matrix));
|
||||
double *elements_data = (double*)malloc(rank*rank*order*sizeof(double));
|
||||
gsl_matrix *inverses = (gsl_matrix*)malloc(order*sizeof(gsl_matrix));
|
||||
double *inverses_data = (double*)malloc(rank*rank*order*sizeof(double));
|
||||
gsl_vector *i1 = gsl_vector_alloc(rank);
|
||||
gsl_vector *i2 = gsl_vector_alloc(rank);
|
||||
gsl_matrix *current_element = gsl_matrix_alloc(rank, rank);
|
||||
queue_t queue;
|
||||
int current;
|
||||
|
||||
|
||||
for(int i = 0; i < rank; i++)
|
||||
generators[i] = gsl_matrix_view_array(generators_data + i*rank*rank, rank, rank).matrix;
|
||||
for(int i = 0; i < order; i++)
|
||||
elements[i] = gsl_matrix_view_array(elements_data + i*rank*rank, rank, rank).matrix;
|
||||
for(int i = 0; i < order; i++)
|
||||
inverses[i] = gsl_matrix_view_array(inverses_data + i*rank*rank, rank, rank).matrix;
|
||||
|
||||
generate_schlaefli_matrix(type, schlaefliMatrix);
|
||||
generate_starting_vector(rank, schlaefliMatrix, startingVector);
|
||||
|
||||
for(int i = 0; i < rank; i++) {
|
||||
gsl_matrix_set_identity(&generators[i]);
|
||||
for(int j = 0; j < rank; j++)
|
||||
if(i == j)
|
||||
gsl_matrix_set(&generators[i], i, j, -1.0);
|
||||
else
|
||||
gsl_matrix_set(&generators[i], i, j, -gsl_matrix_get(schlaefliMatrix, i, j));
|
||||
// gsl_matrix_fprintf(stdout, &generators[i], "%f");
|
||||
// printf("\n");
|
||||
}
|
||||
|
||||
queue_init(&queue);
|
||||
queue_put(&queue, 0);
|
||||
element_count = 1;
|
||||
gsl_matrix_set_identity(&elements[0]);
|
||||
gsl_matrix_set_identity(&inverses[0]);
|
||||
while((current = queue_get(&queue)) != -1) {
|
||||
for(int i = 0; i < rank; i++) {
|
||||
int j;
|
||||
for(j = 0; j < element_count; j++) {
|
||||
if(check_equal_sector(&generators[i], &elements[current], &inverses[j], schlaefliMatrix, startingVector, i1, i2)) { // generators[i] * elements[current] = elements[j]
|
||||
result[rank*current + i] = j;
|
||||
break;
|
||||
}
|
||||
}
|
||||
// if no existing element equals generators[i] * elements[current], create one
|
||||
if(j == element_count) {
|
||||
ERROR(element_count >= order, "Got more elements than the order of the group should be!\n");
|
||||
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &generators[i], &elements[current], 0.0, &elements[element_count]);
|
||||
invert(&elements[element_count], &inverses[element_count]);
|
||||
result[rank*current + i] = element_count;
|
||||
queue_put(&queue, element_count);
|
||||
element_count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
ERROR(element_count != order, "Somethink went wrong building the Coxeter group. Found %d elements, %d expected\n", element_count, order);
|
||||
|
||||
/*
|
||||
for(int i = 0; i < order; i++) {
|
||||
printf("%d: ", i);
|
||||
for(int j = 0; j < rank; j++)
|
||||
printf("%d ", result[rank*i + j]);
|
||||
printf("\n");
|
||||
}
|
||||
*/
|
||||
|
||||
gsl_vector_free(startingVector);
|
||||
gsl_vector_free(i1);
|
||||
gsl_vector_free(i2);
|
||||
gsl_matrix_free(schlaefliMatrix);
|
||||
gsl_matrix_free(current_element);
|
||||
free(generators);
|
||||
free(elements);
|
||||
free(generators_data);
|
||||
free(elements_data);
|
||||
}
|
||||
18
coxeter.h
Normal file
18
coxeter.h
Normal file
@ -0,0 +1,18 @@
|
||||
#ifndef COXETER_H
|
||||
#define COXETER_H
|
||||
|
||||
typedef struct {
|
||||
char series;
|
||||
int rank;
|
||||
} simple_type_t;
|
||||
|
||||
typedef struct {
|
||||
int n;
|
||||
simple_type_t *factors;
|
||||
} semisimple_type_t;
|
||||
|
||||
void generate_coxeter_graph(semisimple_type_t type, int *result);
|
||||
int coxeter_order(semisimple_type_t type);
|
||||
int coxeter_rank(semisimple_type_t type);
|
||||
|
||||
#endif
|
||||
41
queue.h
Normal file
41
queue.h
Normal file
@ -0,0 +1,41 @@
|
||||
#ifndef QUEUE_H
|
||||
#define QUEUE_H
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#define QUEUE_SIZE 2000
|
||||
|
||||
#define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);}
|
||||
|
||||
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
|
||||
467
thickenings.c
Normal file
467
thickenings.c
Normal file
@ -0,0 +1,467 @@
|
||||
#define _GNU_SOURCE
|
||||
|
||||
#include <stdio.h>
|
||||
#include <limits.h>
|
||||
#include <stdlib.h>
|
||||
#include <malloc.h>
|
||||
#include <memory.h>
|
||||
|
||||
#include "coxeter.h"
|
||||
#include "queue.h"
|
||||
|
||||
#define DEBUG(msg, ...) do{fprintf(stderr, msg, ##__VA_ARGS__); }while(0)
|
||||
|
||||
typedef struct _edgelist {
|
||||
int to;
|
||||
struct _edgelist *next;
|
||||
} edgelist_t;
|
||||
|
||||
typedef struct {
|
||||
int *word;
|
||||
int wordlength;
|
||||
int *left;
|
||||
int *right;
|
||||
int opposite;
|
||||
edgelist_t *bruhat_lower;
|
||||
edgelist_t *bruhat_higher;
|
||||
int is_hyperplane_reflection; // boolean value
|
||||
} node_t;
|
||||
|
||||
static char *alphabetize(int *word, int len, const char *alphabet, char *buffer)
|
||||
{
|
||||
int i = 0;
|
||||
for(i = 0; i < len; i++)
|
||||
buffer[i] = alphabet[word[i]];
|
||||
buffer[i] = 0;
|
||||
|
||||
return buffer;
|
||||
}
|
||||
|
||||
static int compare_wordlength(const void *a, const void *b, void *gr)
|
||||
{
|
||||
int i = *((int*)a);
|
||||
int j = *((int*)b);
|
||||
node_t *graph = (node_t*)gr;
|
||||
|
||||
return graph[i].wordlength - graph[j].wordlength;
|
||||
}
|
||||
|
||||
int main(int argc, const char *argv[])
|
||||
{
|
||||
queue_t queue;
|
||||
|
||||
// heap stuff
|
||||
node_t *graph, *graph_unsorted;
|
||||
int *graph_data;
|
||||
int *wordlength_order, *reverse_wordlength_order, *seen, *level;
|
||||
int *words;
|
||||
edgelist_t *edgelists;
|
||||
int *left, *right;
|
||||
int *left_invariant, *right_invariant;
|
||||
|
||||
edgelist_t *edge, *previous;
|
||||
int rank, order;
|
||||
semisimple_type_t type;
|
||||
int edgelist_count, hyperplane_count, max_wordlength;
|
||||
int current, head, i, current_level;
|
||||
int is_fat, is_slim;
|
||||
int thickenings_count, fat_count, slim_count, balanced_count;
|
||||
|
||||
char *string_buffer1, *string_buffer2;
|
||||
|
||||
const char *alphabet = "abcdefghijklmnopqrstuvwxyz";
|
||||
|
||||
ERROR(argc < 2, "Too few arguments!\n");
|
||||
|
||||
type.n = argc - 1;
|
||||
type.factors = (simple_type_t*)malloc((argc-1)*sizeof(simple_type_t));
|
||||
for(int i = 0; i < argc - 1; i++) {
|
||||
type.factors[i].series = argv[i+1][0];
|
||||
type.factors[i].rank = argv[i+1][1] - '0';
|
||||
ERROR(argv[i+1][0] < 'A' || argv[i+1][0] > 'I' || argv[i+1][1] < '1' || argv[i+1][1] > '9', "Arguments must be Xn with X out of A-I and n out of 0-9\n");
|
||||
}
|
||||
|
||||
rank = coxeter_rank(type);
|
||||
order = coxeter_order(type);
|
||||
|
||||
ERROR(strlen(alphabet) < rank, "The alphabet has too few letters\n");
|
||||
|
||||
DEBUG("The group has rank %d and order %d\n", rank, order);
|
||||
|
||||
graph = (node_t*)malloc(order*sizeof(node_t));
|
||||
graph_unsorted = (node_t*)malloc(order*sizeof(node_t));
|
||||
graph_data = (int*)malloc(order*rank*sizeof(int));
|
||||
wordlength_order = (int*)malloc(order*sizeof(int));
|
||||
reverse_wordlength_order = (int*)malloc(order*sizeof(int));
|
||||
seen = (int*)malloc(order*sizeof(int));
|
||||
level = (int*)malloc(order*sizeof(int));
|
||||
left = (int*)malloc(order*rank*sizeof(int));
|
||||
right = (int*)malloc(order*rank*sizeof(int));
|
||||
left_invariant = (int*)malloc(rank*sizeof(int));
|
||||
right_invariant = (int*)malloc(rank*sizeof(int));
|
||||
|
||||
DEBUG("Generate Cayley graph\n");
|
||||
|
||||
generate_coxeter_graph(type, graph_data);
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
graph_unsorted[i].left = &left[i*rank];
|
||||
graph_unsorted[i].right = &right[i*rank];
|
||||
for(int j = 0; j < rank; j++)
|
||||
graph_unsorted[i].left[j] = graph_data[i*rank + j];
|
||||
graph_unsorted[i].word = 0;
|
||||
graph_unsorted[i].wordlength = INT_MAX;
|
||||
graph_unsorted[i].bruhat_lower = 0;
|
||||
graph_unsorted[i].bruhat_higher = 0;
|
||||
graph_unsorted[i].is_hyperplane_reflection = 0;
|
||||
}
|
||||
|
||||
DEBUG("Find wordlengths\n");
|
||||
|
||||
graph_unsorted[0].wordlength = 0;
|
||||
queue_init(&queue);
|
||||
queue_put(&queue, 0);
|
||||
while((current = queue_get(&queue)) != -1) {
|
||||
for(int i = 0; i < rank; i++) {
|
||||
int neighbor = graph_unsorted[current].left[i];
|
||||
if(graph_unsorted[neighbor].wordlength > graph_unsorted[current].wordlength + 1) {
|
||||
graph_unsorted[neighbor].wordlength = graph_unsorted[current].wordlength + 1;
|
||||
queue_put(&queue, neighbor);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
max_wordlength = 0;
|
||||
for(int i = 0; i < order; i++)
|
||||
if(graph_unsorted[i].wordlength > max_wordlength)
|
||||
max_wordlength = graph_unsorted[i].wordlength;
|
||||
|
||||
string_buffer1 = (char*)malloc((max_wordlength+1)*sizeof(char));
|
||||
string_buffer2 = (char*)malloc((max_wordlength+1)*sizeof(char));
|
||||
|
||||
DEBUG("Sort by wordlength\n");
|
||||
|
||||
for(int i = 0; i < order; i++)
|
||||
wordlength_order[i] = i;
|
||||
qsort_r(wordlength_order, order, sizeof(int), compare_wordlength, graph_unsorted); // so wordlength_order is a map new index -> old index
|
||||
for(int i = 0; i < order; i++)
|
||||
reverse_wordlength_order[wordlength_order[i]] = i; // reverse_wordlength_order is a map old index -> new index
|
||||
for(int i = 0; i < order; i++) {
|
||||
graph[i] = graph_unsorted[wordlength_order[i]]; // copy the whole thing
|
||||
for(int j = 0; j < rank; j++)
|
||||
graph[i].left[j] = reverse_wordlength_order[graph[i].left[j]]; // rewrite references
|
||||
}
|
||||
|
||||
DEBUG("Find words\n");
|
||||
|
||||
words = (int*)malloc(order*max_wordlength*sizeof(int));
|
||||
memset(words, 0, order*max_wordlength*sizeof(int));
|
||||
graph[0].word = &words[0];
|
||||
queue_init(&queue);
|
||||
queue_put(&queue, 0);
|
||||
while((current = queue_get(&queue)) != -1) {
|
||||
for(int i = 0; i < rank; i++) {
|
||||
int neighbor = graph[current].left[i];
|
||||
if(graph[neighbor].wordlength == graph[current].wordlength + 1 && graph[neighbor].word == 0) {
|
||||
graph[neighbor].word = &words[neighbor*max_wordlength];
|
||||
memcpy(&graph[neighbor].word[1], &graph[current].word[0], graph[current].wordlength*sizeof(int));
|
||||
graph[neighbor].word[0] = i;
|
||||
queue_put(&queue, neighbor);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
DEBUG("Generate right edges\n");
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
for(int j = 0; j < rank; j++) {
|
||||
current = graph[0].left[j];
|
||||
for(int k = graph[i].wordlength - 1; k >= 0; k--) { // apply group element from right to left
|
||||
current = graph[current].left[graph[i].word[k]];
|
||||
}
|
||||
graph[i].right[j] = current;
|
||||
}
|
||||
}
|
||||
|
||||
DEBUG("Find opposites\n");
|
||||
|
||||
node_t *longest = &graph[order-1];
|
||||
for(int i = 0; i < order; i++) {
|
||||
current = i;
|
||||
for(int k = longest->wordlength - 1; k >= 0; k--)
|
||||
current = graph[current].left[longest->word[k]];
|
||||
graph[i].opposite = current;
|
||||
}
|
||||
|
||||
DEBUG("Enumerate hyperplanes\n"); // every right edge is a reflection along a hyperplane; calculate what this reflection does to the identity
|
||||
|
||||
hyperplane_count = 0;
|
||||
for(int i = 0; i < order; i++) {
|
||||
for(int j = 0; j < rank; j++) {
|
||||
current = 0;
|
||||
int *word1 = graph[i].word;
|
||||
int word1len = graph[i].wordlength;
|
||||
int *word2 = graph[graph[i].right[j]].word; // want to calculate word2 * word1^{-1}
|
||||
int word2len = graph[graph[i].right[j]].wordlength;
|
||||
for(int k = 0; k < word1len; k++) // apply inverse, i.e. go from left to right
|
||||
current = graph[current].left[word1[k]];
|
||||
for(int k = word2len - 1; k >= 0; k--) // now from right to left
|
||||
current = graph[current].left[word2[k]];
|
||||
if(graph[current].is_hyperplane_reflection == 0) {
|
||||
graph[current].is_hyperplane_reflection = 1;
|
||||
hyperplane_count++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
DEBUG("We have %d hyperplanes\n", hyperplane_count);
|
||||
DEBUG("Generate folding order\n");
|
||||
|
||||
edgelists = (edgelist_t*)malloc(order*hyperplane_count*sizeof(edgelist_t));
|
||||
for(int i = 0; i < order; i++) {
|
||||
if(graph[i].is_hyperplane_reflection) {
|
||||
for(int j = 0; j < order; j++) {
|
||||
|
||||
current = j;
|
||||
for(int k = graph[i].wordlength - 1; k >= 0; k--) // apply hyperplane reflection
|
||||
current = graph[current].left[graph[i].word[k]];
|
||||
|
||||
if(graph[j].wordlength < graph[current].wordlength) { // current has higher bruhat order than j
|
||||
edgelists[edgelist_count].to = j;
|
||||
edgelists[edgelist_count].next = graph[current].bruhat_lower;
|
||||
graph[current].bruhat_lower = &edgelists[edgelist_count];
|
||||
edgelist_count++;
|
||||
} else if(graph[j].wordlength > graph[current].wordlength) { // j has higher bruhat order than current; these are already included from the other side
|
||||
} else {
|
||||
ERROR(1, "Chambers of equal word lengths should not be folded on each other!\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
DEBUG("Remove redundant edges\n");
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
memset(seen, 0, order*sizeof(int));
|
||||
for(int len = 1; len <= max_wordlength; len++) {
|
||||
// remove all edges originating from i of length len which connect to something already seen using shorter edges
|
||||
edge = graph[i].bruhat_lower;
|
||||
previous = (edgelist_t*)0;
|
||||
while(edge) {
|
||||
if(seen[edge->to] && graph[i].wordlength - graph[edge->to].wordlength == len) {
|
||||
// printf("deleting from %d to %d\n", i, edge->to);
|
||||
if(previous)
|
||||
previous->next = edge->next;
|
||||
else
|
||||
graph[i].bruhat_lower = edge->next;
|
||||
} else {
|
||||
previous = edge;
|
||||
}
|
||||
edge = edge->next;
|
||||
}
|
||||
|
||||
// see which nodes we can reach using only edges up to length len, mark them as seen
|
||||
queue_init(&queue);
|
||||
queue_put(&queue, i);
|
||||
seen[i] = 1;
|
||||
while((current = queue_get(&queue)) != -1) {
|
||||
edge = graph[current].bruhat_lower;
|
||||
while(edge) {
|
||||
if(!seen[edge->to] && graph[current].wordlength - graph[edge->to].wordlength == len) {
|
||||
seen[edge->to] = 1;
|
||||
queue_put(&queue, edge->to);
|
||||
}
|
||||
edge = edge->next;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
DEBUG("Reverse folding order\n");
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
edge = graph[i].bruhat_lower;
|
||||
while(edge) {
|
||||
edgelists[edgelist_count].to = i;
|
||||
edgelists[edgelist_count].next = graph[edge->to].bruhat_higher;
|
||||
graph[edge->to].bruhat_higher = &edgelists[edgelist_count];
|
||||
edgelist_count++;
|
||||
edge = edge->next;
|
||||
}
|
||||
}
|
||||
|
||||
printf("The group is: \n");
|
||||
for(int i =0; i < order; i++) {
|
||||
if(graph[i].wordlength == 0)
|
||||
printf("1 ");
|
||||
else
|
||||
printf("%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1));
|
||||
}
|
||||
printf("\n");
|
||||
|
||||
DEBUG("Enumerate thickenings\n");
|
||||
DEBUG("\n");
|
||||
|
||||
thickenings_count = fat_count = slim_count = balanced_count = 0;
|
||||
memset(level, 0, order*sizeof(int));
|
||||
current_level = 1;
|
||||
head = order - 1;
|
||||
level[head] = -1;
|
||||
while(current_level > 0) {
|
||||
// calculate transitive closure
|
||||
queue_init(&queue);
|
||||
queue_put(&queue, head);
|
||||
while((current = queue_get(&queue)) != -1) {
|
||||
edge = graph[current].bruhat_lower;
|
||||
while(edge) {
|
||||
if(level[edge->to] == 0) {
|
||||
level[edge->to] = current_level;
|
||||
queue_put(&queue, edge->to);
|
||||
}
|
||||
edge = edge->next;
|
||||
}
|
||||
}
|
||||
|
||||
is_fat = is_slim = 1;
|
||||
for(int i = 0; i < order; i++) {
|
||||
if(level[graph[i].opposite] != 0) {
|
||||
if(level[i] != 0)
|
||||
is_slim = 0;
|
||||
} else {
|
||||
if(level[i] == 0)
|
||||
is_fat = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// count
|
||||
thickenings_count++;
|
||||
if(is_fat)
|
||||
fat_count++;
|
||||
if(is_slim)
|
||||
slim_count++;
|
||||
if(is_slim && is_fat)
|
||||
balanced_count++;
|
||||
|
||||
// print out levels
|
||||
if(is_fat && is_slim) {
|
||||
// check for invariances
|
||||
for(int j = 0; j < rank; j++) {
|
||||
left_invariant[j] = 1;
|
||||
right_invariant[j] = 1;
|
||||
}
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
for(int j = 0; j < rank; j++) {
|
||||
if(level[i] == 0 && level[graph[i].left[j]] != 0 || level[i] != 0 && level[graph[i].left[j]] == 0)
|
||||
left_invariant[j] = 0;
|
||||
if(level[i] == 0 && level[graph[i].right[j]] != 0 || level[i] != 0 && level[graph[i].right[j]] == 0)
|
||||
right_invariant[j] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
if(level[i] != 0)
|
||||
printf("x");
|
||||
else
|
||||
printf("%d", level[i]);
|
||||
}
|
||||
printf(" left: ");
|
||||
for(int j = 0; j < rank; j++)
|
||||
if(left_invariant[j])
|
||||
printf("%c", alphabet[j]);
|
||||
else
|
||||
printf(" ");
|
||||
printf(" right: ");
|
||||
for(int j = 0; j < rank; j++)
|
||||
if(right_invariant[j])
|
||||
printf("%c", alphabet[j]);
|
||||
else
|
||||
printf(" ");
|
||||
printf("\n");
|
||||
}
|
||||
|
||||
// try to find empty spot to the left of "head"
|
||||
for(i = head - 1; i >= 0; i--)
|
||||
if(level[i] == 0)
|
||||
break;
|
||||
if(i >= 0) {
|
||||
head = i;
|
||||
level[head] = -1;
|
||||
current_level++;
|
||||
continue;
|
||||
}
|
||||
|
||||
// if none was found, try to move "head" to the left
|
||||
while(current_level > 0) {
|
||||
for(i = head - 1; i >= 0; i--)
|
||||
if(level[i] == 0 || level[i] >= current_level)
|
||||
break;
|
||||
if(i >= 0) { // if this was successful, just move head
|
||||
level[head] = 0;
|
||||
head = i;
|
||||
level[head] = -1;
|
||||
break;
|
||||
} else { // if moving the head is not possible, take the next head to the right
|
||||
current_level--;
|
||||
level[head] = 0;
|
||||
do {
|
||||
head++;
|
||||
} while(head < order && level[head] != -1);
|
||||
}
|
||||
}
|
||||
|
||||
// clean up
|
||||
for(int i = 0; i < head; i++)
|
||||
if(level[i] >= current_level)
|
||||
level[i] = 0;
|
||||
|
||||
// print out levels
|
||||
/* for(int i = 0; i < order; i++) { */
|
||||
/* if(level[i] == -1) */
|
||||
/* printf("x"); */
|
||||
/* else */
|
||||
/* printf("%d", level[i]); */
|
||||
/* } */
|
||||
/* printf(" %d\n", current_level); */
|
||||
}
|
||||
|
||||
printf("\n");
|
||||
printf("Found %d thickenings, %d fat, %d slim, %d balanced\n", thickenings_count, fat_count, slim_count, balanced_count);
|
||||
|
||||
/*
|
||||
printf("\n");
|
||||
|
||||
for(int i = 0; i < order; i++) {
|
||||
printf("%-6s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1));
|
||||
|
||||
edge = graph[i].bruhat_higher;
|
||||
while(edge) {
|
||||
printf("-> %-6s ", alphabetize(graph[edge->to].word, graph[edge->to].wordlength, alphabet, string_buffer1));
|
||||
edge = edge->next;
|
||||
}
|
||||
//printf("Bruhat: %6s -> %-6s\n", alphabetize(graph[current].word, graph[current].wordlength, alphabet, string_buffer1), alphabetize(graph[j].word, graph[j].wordlength, alphabet, string_buffer2));
|
||||
|
||||
printf("\n");
|
||||
}
|
||||
printf("\n");
|
||||
*/
|
||||
|
||||
free(edgelists);
|
||||
free(words);
|
||||
free(string_buffer1);
|
||||
free(string_buffer2);
|
||||
free(graph);
|
||||
free(graph_unsorted);
|
||||
free(graph_data);
|
||||
free(wordlength_order);
|
||||
free(reverse_wordlength_order);
|
||||
free(seen);
|
||||
free(level);
|
||||
free(left);
|
||||
free(right);
|
||||
free(left_invariant);
|
||||
free(right_invariant);
|
||||
free(type.factors);
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user