diff --git a/Makefile b/Makefile index 7b6c37c..b3eafce 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,6 @@ HEADERS=coxeter.h thickenings.h queue.h OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99 +#OPTIONS=-O0 -g -std=gnu99 all: generate process diff --git a/coxeter.c b/coxeter.c index 10a9e96..ca8ab50 100644 --- a/coxeter.c +++ b/coxeter.c @@ -176,6 +176,37 @@ int coxeter_hyperplanes(semisimple_type_t type) return hyperplanes; } +unsigned long opposition_involution(semisimple_type_t type, unsigned long theta) +{ + int offset = 0; + unsigned long result = 0; + for(int i = 0; i < type.n; i++) { + unsigned long current = (theta >> offset) & ((1 << type.factors[i].rank) - 1); + unsigned long iota_current; + + if(type.factors[i].series == 'B' || type.factors[i].series == 'C' || type.factors[i].series == 'F' || type.factors[i].series == 'G') { + iota_current = current; + } else if(type.factors[i].series == 'A') { + iota_current = 0; + for(int j = 0; j < type.factors[i].rank; j++) + iota_current += ((current >> j) & 1) << (type.factors[i].rank - 1 - j); + } else if(type.factors[i].series == 'D') { + if(type.factors[i].rank % 2 == 0) { + iota_current = current; + } else { + ERROR(1, "The opposition involution for type %c%d is not yet implemented!\n", type.factors[i].series, type.factors[i].rank); + } + } else if(type.factors[i].series == 'E') { + ERROR(1, "The opposition involution for En is not yet implemented!\n"); + } + + result += iota_current << offset; + offset += type.factors[i].rank; + } + + return result; +} + static void generate_starting_vector(int rank, gsl_matrix *schlaefli, gsl_vector *result) { gsl_matrix *schlaefliLU = gsl_matrix_alloc(rank, rank); diff --git a/coxeter.h b/coxeter.h index babe2e9..b62278e 100644 --- a/coxeter.h +++ b/coxeter.h @@ -15,5 +15,6 @@ void generate_coxeter_graph(semisimple_type_t type, int *result); int coxeter_order(semisimple_type_t type); int coxeter_hyperplanes(semisimple_type_t type); int coxeter_rank(semisimple_type_t type); +unsigned long opposition_involution(semisimple_type_t type, unsigned long theta); #endif diff --git a/generate.c b/generate.c index 6e4c87f..818ac16 100644 --- a/generate.c +++ b/generate.c @@ -8,56 +8,62 @@ int main(int argc, const char *argv[]) { semisimple_type_t type; + unsigned long right_invariance, left_invariance; + int rank, order, hyperplanes, cosets; - // heap stuff node_t *graph; - int *left, *right; - edgelist_t *edgelists; - int *words; - - int rank, order; char string_buffer1[1000]; const char *alphabet = "abcdefghijklmnopqrstuvwxyz"; + // read arguments + ERROR(argc < 2, "Too few arguments!\n"); - type.n = argc - 1; - type.factors = (simple_type_t*)malloc((argc-1)*sizeof(simple_type_t)); + type.n = 0; for(int i = 0; i < argc - 1; i++) { + if(argv[i+1][0] < 'A' || argv[i+1][0] > 'I') + break; + type.n++; + } + + type.factors = (simple_type_t*)malloc(type.n*sizeof(simple_type_t)); + for(int i = 0; i < type.n; 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); + left_invariance = right_invariance = 0; + + if(argc - type.n >= 3) { + if(strcmp(argv[type.n + 1], "-") != 0) + for(int i = 0; i < strlen(argv[type.n + 1]); i++) + left_invariance |= (1 << (argv[type.n + 1][i] - 'a')); + if(strcmp(argv[type.n + 2], "-") != 0) + for(int i = 0; i < strlen(argv[type.n + 2]); i++) + right_invariance |= (1 << (argv[type.n + 2][i] - 'a')); + } ERROR(strlen(alphabet) < rank, "The alphabet has too few letters\n"); - // initialize - - graph = (node_t*)malloc(order*sizeof(node_t)); - left = (int*)malloc(order*rank*sizeof(int)); - right = (int*)malloc(order*rank*sizeof(int)); - - for(int i = 0; i < order; i++) { - graph[i].left = &left[rank*i]; - graph[i].right = &right[rank*i]; - } - // generate graph - prepare_graph(type, graph, &edgelists, &words); + graph = graph_alloc(type); + cosets = prepare_simplified_graph(type, left_invariance, right_invariance, graph); + + ERROR(cosets < 0, "The left invariance is not preserved by the opposition involution: %d %d!\n", left_invariance, opposition_involution(type, left_invariance)); // print stuff - int hyperplane_count = coxeter_hyperplanes(type); + rank = coxeter_rank(type); // number of simple roots + order = coxeter_order(type); // number of Weyl group elements + hyperplanes = coxeter_hyperplanes(type); // number of positive roots - fprintf(stderr, "Rank: %d\t\tOrder: %d\t\tHyperplanes: %d\n", rank, order, hyperplane_count); + fprintf(stderr, "Rank: %d\tOrder: %d\tPositive Roots: %d\tCosets: %d\n", rank, order, hyperplanes, cosets); fprintf(stderr, "\n"); - fprintf(stderr, "Group elements: \n"); - for(int i = 0, wl = 0; i < order; i++) { + fprintf(stderr, "Shortest coset representatives: \n"); + for(int i = 0, wl = 0; i < cosets; i++) { if(i == 0) { fprintf(stderr, "1"); } else if(graph[i].wordlength > wl) { @@ -72,16 +78,15 @@ int main(int argc, const char *argv[]) fwrite(&type.n, sizeof(int), 1, stdout); fwrite(type.factors, sizeof(simple_type_t), type.n, stdout); - long count = enumerate_balanced_thickenings(type, graph, order, alphabet, stdout); + fwrite(&left_invariance, sizeof(unsigned long), type.n, stdout); + fwrite(&right_invariance, sizeof(unsigned long), type.n, stdout); + + long count = enumerate_balanced_thickenings(type, graph, cosets, alphabet, stdout); fprintf(stderr, "\n"); fprintf(stderr, "Found %ld balanced thickenings\n\n", count); - free(graph); - free(left); - free(right); - free(edgelists); - free(words); + graph_free(type, graph); free(type.factors); return 0; diff --git a/process-old.c b/process-old.c new file mode 100644 index 0000000..0dd3fb4 --- /dev/null +++ b/process-old.c @@ -0,0 +1,193 @@ +#include +#include +#include + +#include "thickenings.h" +#include "coxeter.h" +#include "queue.h" + +#define SWAP(t, a, b) do {t tmp = a; a = b; b = tmp;} while(0) + +char *stringify_SLn1_permutation(int *word, int wordlength, int rank, char *str) +{ + for(int i = 0; i <= rank; i++) + str[i] = '1' + i; + str[rank+1] = 0; + + for(int i = 0; i < wordlength; i++) { + char tmp = str[word[i]]; + str[word[i]] = str[word[i]+1]; + str[word[i]+1] = tmp; + } + return str; +} + + +char *stringify_Onn1_permutation(int *word, int wordlength, int rank, char *str) +{ + for(int i = 0; i <= rank*2; i++) + str[i] = '1' + i; + str[2*rank+1] = 0; + + for(int i = 0; i < wordlength; i++) { + if(word[i] == 0) + SWAP(char, str[rank-1], str[rank+1]); + else { + SWAP(char, str[rank-word[i]], str[rank-word[i]-1]); + SWAP(char, str[rank+word[i]], str[rank+word[i]+1]); + } + } + return str; +} + +int main(int argc, const char *argv[]) +{ + FILE *infile; + struct stat st; + int rank, order, hyperplanes; + semisimple_type_t type; + int n; + signed char *level; + node_t *graph; + int *left, *right; + int left_invariant, right_invariant; + int left_invariant_wanted = 0, right_invariant_wanted = 0; + + unsigned long left_invariance, right_invariance; + edgelist_t *edgelists; + int *words; + + queue_t queue; + int current; + int *seen; + int *generators; + int ngens; + + char string_buffer1[1000]; + const char *alphabet = "abcdefghijklmnopqrstuvwxyz"; + + // parse arguments + + if(argc < 2) + infile = stdin; + else { + if(strcmp(argv[1], "-") == 0) + infile = stdin; + else + infile = fopen(argv[1], "rb"); + + if(argc >= 4) { + if(strcmp(argv[2], "-") != 0) + for(int i = 0; i < strlen(argv[2]); i++) + left_invariant_wanted |= (1 << (argv[2][i] - 'a')); + if(strcmp(argv[3], "-") != 0) + for(int i = 0; i < strlen(argv[3]); i++) + right_invariant_wanted |= (1 << (argv[3][i] - 'a')); + } + } + + fread(&type.n, sizeof(int), 1, infile); // we completely trust the input data + type.factors = malloc(type.n * sizeof(simple_type_t)); + fread(type.factors, sizeof(simple_type_t), type.n, infile); + fread(&left_invariance, sizeof(simple_type_t), type.n, infile); + fread(&right_invariance, sizeof(simple_type_t), type.n, infile); + + // get graph + + rank = coxeter_rank(type); + order = coxeter_order(type); + hyperplanes = coxeter_hyperplanes(type); + ERROR(strlen(alphabet) < rank, "The alphabet has too few letters\n"); + seen = (int*)malloc(order*sizeof(int)); + generators = (int*)malloc(order*sizeof(int)); + level = (signed char*)malloc(order*sizeof(int)); + + graph = graph_alloc(type); + prepare_graph(type, graph); + + // finally do stuff + + int counter = 0; + + while(fread(level, sizeof(signed char), order, infile) == order) { + + /* + if((counter++) % 100000 == 0) + print_thickening(rank, order, level, 0, 0, 0, alphabet, stdout); + continue; + */ + + left_invariant = right_invariant = -1; // all 1s + for(int j = 0; j < order; j++) { + for(int k = 0; k < rank; k++) { + if(level[j] > 0 && level[graph[j].left[k]] < 0 || level[j] < 0 && level[graph[j].left[k]] > 0) { + left_invariant &= ~(1 << k); + } + if(level[j] > 0 && level[graph[j].right[k]] < 0 || level[j] < 0 && level[graph[j].right[k]] > 0) { + right_invariant &= ~(1 << k); + } + } + } + + if((~left_invariant & left_invariant_wanted) == 0 && (~right_invariant & right_invariant_wanted) == 0) { + ngens = 0; + memset(generators, 0, order*sizeof(int)); + for(int j = 0; j < order; j++) { + if(level[j] == HEAD_MARKER && generators[j] == 0) { // ignore the generator, if it is equivalent to one already seen + ngens++; + queue_init(&queue); + queue_put(&queue, j); + while((current = queue_get(&queue)) != -1) { + if(generators[current] == 0) { // visit everyone only once + generators[current] = ngens; + for(int k = 0; k < rank; k++) { + if(left_invariant & (1 << k)) + queue_put(&queue, graph[current].left[k]); + if(right_invariant & (1 << k)) + queue_put(&queue, graph[current].right[k]); + } + } + } + } + } + + printf("left: "); + for(int j = 0; j < rank; j++) + printf("%c", left_invariant & (1 << j) ? alphabet[j] : ' '); + printf(" right: "); + for(int j = 0; j < rank; j++) + printf("%c", right_invariant & (1 << j) ? alphabet[j] : ' '); + printf(" generators: "); + + memset(seen, 0, order*sizeof(int)); + for(int i = 0; i < order; i++) { + if(generators[i] != 0 && seen[generators[i]-1] == 0) { + seen[generators[i]-1] = 1; + // if(type.n == 1 && type.factors[0].series == 'A') + // printf("%s ", stringify_SLn1_permutation(graph[i].word, graph[i].wordlength, rank, string_buffer1)); + // else if(type.n == 1 && (type.factors[0].series == 'B' || type.factors[0].series == 'C')) + // printf("%s ", stringify_Onn1_permutation(graph[i].word, graph[i].wordlength, rank, string_buffer1)); + // else + if(i == 0) + printf("1 "); + else + printf("%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + } + } + + printf("\n"); + } + } + + if(infile != stdin) + fclose(infile); + + // cleanup + + graph_free(type, graph); + free(seen); + free(generators); + free(type.factors); + + return 0; +} diff --git a/process.c b/process.c index c62cad2..0c9759a 100644 --- a/process.c +++ b/process.c @@ -6,105 +6,117 @@ #include "coxeter.h" #include "queue.h" -#define SWAP(t, a, b) do {t tmp = a; a = b; b = tmp;} while(0) - -char *stringify_SLn1_permutation(int *word, int wordlength, int rank, char *str) -{ - for(int i = 0; i <= rank; i++) - str[i] = '1' + i; - str[rank+1] = 0; - - for(int i = 0; i < wordlength; i++) { - char tmp = str[word[i]]; - str[word[i]] = str[word[i]+1]; - str[word[i]+1] = tmp; - } - return str; -} - - -char *stringify_Onn1_permutation(int *word, int wordlength, int rank, char *str) -{ - for(int i = 0; i <= rank*2; i++) - str[i] = '1' + i; - str[2*rank+1] = 0; - - for(int i = 0; i < wordlength; i++) { - if(word[i] == 0) - SWAP(char, str[rank-1], str[rank+1]); - else { - SWAP(char, str[rank-word[i]], str[rank-word[i]-1]); - SWAP(char, str[rank+word[i]], str[rank+word[i]+1]); - } - } - return str; -} - int main(int argc, const char *argv[]) { FILE *infile; - struct stat st; - int rank, order; semisimple_type_t type; - int n; - int *thickenings; - signed char *level; + unsigned long left_invariance, right_invariance; // these are the invariances we have already modded out + unsigned long left_invariant, right_invariant; // these are the invariances of the thickening under consideration + int rank, cosets; node_t *graph; - int *left, *right; - int left_invariant, right_invariant; - int left_invariant_wanted = 0, right_invariant_wanted = 0; - edgelist_t *edgelists; - int *words; - + signed char *thickening; + int *seen, *generators; queue_t queue; + int ngenerators; int current; - int *seen; - int *generators; - int ngens; char string_buffer1[1000]; const char *alphabet = "abcdefghijklmnopqrstuvwxyz"; - // parse arguments - if(argc < 2) infile = stdin; - else { - if(strcmp(argv[1], "-") == 0) - infile = stdin; - else - infile = fopen(argv[1], "rb"); + else + infile = fopen(argv[1], "rb"); - if(argc >= 4) { - if(strcmp(argv[2], "-") != 0) - for(int i = 0; i < strlen(argv[2]); i++) - left_invariant_wanted |= (1 << (argv[2][i] - 'a')); - if(strcmp(argv[3], "-") != 0) - for(int i = 0; i < strlen(argv[3]); i++) - right_invariant_wanted |= (1 << (argv[3][i] - 'a')); - } - } - - fread(&type.n, sizeof(int), 1, infile); // we completely trust the input data + // we completely trust the input data + ERROR(fread(&type.n, sizeof(int), 1, infile) == 0, "The input file seems to be empty!\n"); type.factors = malloc(type.n * sizeof(simple_type_t)); fread(type.factors, sizeof(simple_type_t), type.n, infile); - - // get graph + fread(&left_invariance, sizeof(simple_type_t), type.n, infile); + fread(&right_invariance, sizeof(simple_type_t), type.n, infile); rank = coxeter_rank(type); - order = coxeter_order(type); - ERROR(strlen(alphabet) < rank, "The alphabet has too few letters\n"); - graph = (node_t*)malloc(order*sizeof(node_t)); - left = (int*)malloc(order*rank*sizeof(int)); - right = (int*)malloc(order*rank*sizeof(int)); + graph = graph_alloc(type); + cosets = prepare_simplified_graph(type, left_invariance, right_invariance, graph); + + thickening = (signed char*)malloc(cosets*sizeof(signed char)); + generators = (int*)malloc(cosets*sizeof(int)); + seen = (int*)malloc(cosets*sizeof(int)); + + while(fread(thickening, sizeof(signed char), cosets, infile) == cosets) { + + // determine invariances of this thickening + left_invariant = right_invariant = -1; // set all bits to 1 + for(int j = 0; j < cosets; j++) { + for(int k = 0; k < rank; k++) { + if(thickening[j] > 0 && thickening[graph[j].left[k]] < 0 || + thickening[j] < 0 && thickening[graph[j].left[k]] > 0) + left_invariant &= ~(1 << k); + if(thickening[j] > 0 && thickening[graph[j].right[k]] < 0 || + thickening[j] < 0 && thickening[graph[j].right[k]] > 0) + right_invariant &= ~(1 << k); + } + } + + // print this stuff + printf("left: "); + for(int j = 0; j < rank; j++) + printf("%c", left_invariant & (1 << j) ? alphabet[j] : ' '); + printf(" right: "); + for(int j = 0; j < rank; j++) + printf("%c", right_invariant & (1 << j) ? alphabet[j] : ' '); + printf(" generators: "); + + // find a minimal set of weyl group elements such that the union of the ideals generated by their cosets wrt the invariances determined above gives the thickening + // in the first step, mark everything which is equivalent to a "head" by a generator id + ngenerators = 0; + memset(generators, 0, cosets*sizeof(int)); + for(int j = 0; j < cosets; j++) { + if(thickening[j] == HEAD_MARKER && generators[j] == 0) { // ignore the generator, if it is equivalent to one already seen + ngenerators++; + queue_init(&queue); + queue_put(&queue, j); + while((current = queue_get(&queue)) != -1) { + if(generators[current] == 0) { // visit everyone only once + generators[current] = ngenerators; + for(int k = 0; k < rank; k++) { + if(left_invariant & (1 << k)) + queue_put(&queue, graph[current].left[k]); + if(right_invariant & (1 << k)) + queue_put(&queue, graph[current].right[k]); + } + } + } + } + } + + // in the second step, go through the list in ascending word length order and print the first appearance of each generator id + memset(seen, 0, cosets*sizeof(int)); + for(int i = 0; i < cosets; i++) { + if(generators[i] != 0 && seen[generators[i]-1] == 0) { + seen[generators[i]-1] = 1; + printf("%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + } + } + + printf("\n"); + } + + if(infile != stdin) + fclose(infile); + + graph_free(type, graph); + free(type.factors); + free(thickening); +} + /******************************************************************************************* + seen = (int*)malloc(order*sizeof(int)); generators = (int*)malloc(order*sizeof(int)); level = (signed char*)malloc(order*sizeof(int)); - for(int i = 0; i < order; i++) { - graph[i].left = &left[i*rank]; - graph[i].right = &right[i*rank]; - } - prepare_graph(type, graph, &edgelists, &words); + + graph = graph_alloc(type); + prepare_graph(type, graph); // finally do stuff @@ -112,24 +124,6 @@ int main(int argc, const char *argv[]) while(fread(level, sizeof(signed char), order, infile) == order) { - /* - if((counter++) % 100000 == 0) - print_thickening(rank, order, level, 0, 0, 0, alphabet, stdout); - continue; - */ - - left_invariant = right_invariant = -1; // all 1s - for(int j = 0; j < order; j++) { - for(int k = 0; k < rank; k++) { - if(level[j] > 0 && level[graph[j].left[k]] < 0 || level[j] < 0 && level[graph[j].left[k]] > 0) { - left_invariant &= ~(1 << k); - } - if(level[j] > 0 && level[graph[j].right[k]] < 0 || level[j] < 0 && level[graph[j].right[k]] > 0) { - right_invariant &= ~(1 << k); - } - } - } - if((~left_invariant & left_invariant_wanted) == 0 && (~right_invariant & right_invariant_wanted) == 0) { ngens = 0; memset(generators, 0, order*sizeof(int)); @@ -185,15 +179,11 @@ int main(int argc, const char *argv[]) // cleanup - free(thickenings); - free(edgelists); - free(words); - free(graph); - free(left); - free(right); + graph_free(type, graph); free(seen); free(generators); free(type.factors); return 0; } +*/ diff --git a/queue.h b/queue.h index aab9539..a19cc31 100644 --- a/queue.h +++ b/queue.h @@ -4,7 +4,7 @@ #include #include -#define QUEUE_SIZE 2000 +#define QUEUE_SIZE 5000 #define ERROR(condition, msg, ...) if(condition){fprintf(stderr, msg, ##__VA_ARGS__); exit(1);} diff --git a/test.c b/test.c index e9163c4..7d86a67 100644 --- a/test.c +++ b/test.c @@ -24,6 +24,60 @@ edgelist_t *edgelist_add(edgelist_t *list, int new, edgelist_t *storage, int *st return new_link; } +int main(int argc, const char *argv[]) +{ + unsigned long left_invariance = atoi(argv[2]); + unsigned long right_invariance = atoi(argv[3]); + semisimple_type_t type; + type.n = 1; + type.factors = (simple_type_t*)malloc(type.n*sizeof(simple_type_t)); + type.factors[0].series = 'B'; + type.factors[0].rank = atoi(argv[1]); + + char alphabet[] = "abcdefghijklmnopqrstuvwxyz"; + char buffer[1024], buffer2[1024]; + + node_t *graph; + int *leftbuf, *rightbuf; + edgelist_t *edgelists; + int *words; + int rank, order, hyperplanes, cosets; + + // initialize + + rank = coxeter_rank(type); + order = coxeter_order(type); + hyperplanes = coxeter_hyperplanes(type); + + graph = (node_t*)malloc(order*sizeof(node_t)); + leftbuf = (int*)malloc(order*rank*sizeof(int)); + rightbuf = (int*)malloc(order*rank*sizeof(int)); + edgelists = (edgelist_t*)malloc(4*order*hyperplanes*sizeof(edgelist_t)); + words = (int*)malloc(order*hyperplanes*sizeof(int)); + + for(int i = 0; i < order; i++) { + graph[i].left = &leftbuf[rank*i]; + graph[i].right = &rightbuf[rank*i]; + } + + // generate graph + + cosets = prepare_simplified_graph(type, left_invariance, right_invariance, graph, edgelists, words); + + // do something + + fprintf(stderr, "There are %d double cosets.\n", cosets); + + // cleanup + + free(leftbuf); + free(rightbuf); + free(graph); + free(edgelists); + free(words); +} + +/* int main(int argc, const char *argv[]) { semisimple_type_t type; @@ -245,3 +299,4 @@ int main(int argc, const char *argv[]) return 0; } +*/ diff --git a/thickenings.c b/thickenings.c index 4855d62..594171a 100644 --- a/thickenings.c +++ b/thickenings.c @@ -57,24 +57,28 @@ static int compare_wordlength(const void *a, const void *b, void *gr) return graph[i].wordlength - graph[j].wordlength; } -void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists_pointer, int **words_pointer) // the edgelists_pointer and words_pointer arguments are just for freeing afterwards +void prepare_graph(semisimple_type_t type, node_t *graph) { queue_t queue; - int rank, order; + edgelist_t *edgelists_lower, *edgelists_higher; + int rank, order, hyperplanes; edgelist_t *edge, *previous; - int edgelist_count, max_wordlength, hyperplane_count; + int edgelist_count, hyperplane_count; int current; int *graph_data; node_t *graph_unsorted; - int *wordlength_order, *reverse_wordlength_order, *seen, *words; - edgelist_t *edgelists; + int *wordlength_order, *reverse_wordlength_order, *seen; // initialize rank = coxeter_rank(type); order = coxeter_order(type); + hyperplanes = coxeter_hyperplanes(type); + + edgelists_higher = graph[0].bruhat_higher; + edgelists_lower = &graph[0].bruhat_higher[order*hyperplanes/2]; graph_data = (int*)malloc(order*rank*sizeof(int)); graph_unsorted = (node_t*)malloc(order*sizeof(node_t)); @@ -83,13 +87,10 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists seen = (int*)malloc(order*sizeof(int)); for(int i = 0; i < order; i++) { - graph_unsorted[i].left = graph[i].left; - graph_unsorted[i].right = graph[i].right; - 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; + graph[i].bruhat_lower = 0; + graph[i].bruhat_higher = 0; + graph[i].is_hyperplane_reflection = 0; } // get coxeter graph @@ -98,7 +99,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists for(int i = 0; i < order; i++) for(int j = 0; j < rank; j++) - graph_unsorted[i].left[j] = graph_data[i*rank + j]; + graph_unsorted[i].left = &graph_data[i*rank]; // find wordlengths @@ -115,11 +116,6 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists } } - max_wordlength = 0; - for(int i = 0; i < order; i++) - if(graph_unsorted[i].wordlength > max_wordlength) - max_wordlength = graph_unsorted[i].wordlength; - // sort by wordlength for(int i = 0; i < order; i++) @@ -128,23 +124,22 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists 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 + // we have only set left and wordlength so far, so just copy these + graph[i].wordlength = graph_unsorted[wordlength_order[i]].wordlength; for(int j = 0; j < rank; j++) - graph[i].left[j] = reverse_wordlength_order[graph[i].left[j]]; // rewrite references + graph[i].left[j] = reverse_wordlength_order[graph_unsorted[wordlength_order[i]].left[j]]; // rewrite references } // find words - words = (int*)malloc(order*max_wordlength*sizeof(int)); - memset(words, 0, order*max_wordlength*sizeof(int)); - graph[0].word = &words[0]; + for(int i = 0; i < order; i++) + memset(graph[i].word, 0, hyperplanes*sizeof(int)); 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]; + if(graph[neighbor].wordlength == graph[current].wordlength + 1 && graph[neighbor].word[0] == 0) { memcpy(&graph[neighbor].word[1], &graph[current].word[0], graph[current].wordlength*sizeof(int)); graph[neighbor].word[0] = i; queue_put(&queue, neighbor); @@ -197,7 +192,6 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists // generate folding order - edgelists = (edgelist_t*)malloc(order*hyperplane_count*sizeof(edgelist_t)); edgelist_count = 0; for(int i = 0; i < order; i++) { if(graph[i].is_hyperplane_reflection) { @@ -208,9 +202,9 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists 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]; + edgelists_lower[edgelist_count].to = j; + edgelists_lower[edgelist_count].next = graph[current].bruhat_lower; + graph[current].bruhat_lower = &edgelists_lower[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 { @@ -263,20 +257,18 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists // reverse folding order + edgelist_count = 0; 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]; + edgelists_higher[edgelist_count].to = i; + edgelists_higher[edgelist_count].next = graph[edge->to].bruhat_higher; + graph[edge->to].bruhat_higher = &edgelists_higher[edgelist_count]; edgelist_count++; edge = edge->next; } } - *edgelists_pointer = edgelists; - *words_pointer = words; - free(graph_data); free(graph_unsorted); free(wordlength_order); @@ -284,6 +276,264 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists free(seen); } +static int edgelist_contains(edgelist_t *list, int x) { + while(list) { + if(list->to == x) + return 1; + list = list->next; + } + return 0; +} + +static edgelist_t *edgelist_add(edgelist_t *list, int new, edgelist_t *storage, int *storage_index) +{ + edgelist_t *new_link = &storage[*storage_index]; + new_link->next = list; + new_link->to = new; + (*storage_index)++; + return new_link; +} + +int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigned long right, node_t *simplified_graph) +{ + node_t *full_graph; + int edgelists_used; + int rank, order, hyperplanes; + int *reduced, *group, *simplified; + int *seen; + int current; + edgelist_t *edge, *previous; + queue_t queue; + int ncosets; + + if(opposition_involution(type, left) != left) + return -1; + + edgelist_t *edgelists_higher = &simplified_graph[0].bruhat_higher[0]; + edgelist_t *edgelists_lower = &simplified_graph[0].bruhat_higher[order*hyperplanes/2]; + + // get full graph + + full_graph = graph_alloc(type); + prepare_graph(type, full_graph); + + // initialize stuff + + rank = coxeter_rank(type); + order = coxeter_order(type); + hyperplanes = coxeter_hyperplanes(type); + + reduced = (int*)malloc(order*sizeof(int)); + group = (int*)malloc(order*sizeof(int)); + simplified = (int*)malloc(order*sizeof(int)); + for(int i = 0; i < order; i++) { + group[i] = -1; + reduced[i] = i; + } + + // step 1: group + for(int i = 0; i < order; i++) { + if(group[i] != -1) + continue; + + queue_init(&queue); + queue_put(&queue, i); + while((current = queue_get(&queue)) != -1) { + if(group[current] != -1) + continue; + group[current] = i; + + for(int j = 0; j < rank; j++) { + if(left & (1 << j)) + queue_put(&queue, full_graph[current].left[j]); + if(right & (1 << j)) + queue_put(&queue, full_graph[current].right[j]); + } + } + } + + // step 2: find minimum + for(int i = 0; i < order; i++) + if(full_graph[i].wordlength < full_graph[reduced[group[i]]].wordlength) + reduced[group[i]] = i; + + // step 3: assign minimum to all + for(int i = 0; i < order; i++) + reduced[i] = reduced[group[i]]; + + // step 4: assign indices to cosets + ncosets = 0; + for(int i = 0; i < order; i++) + if(reduced[i] == i) + simplified[i] = ncosets++; + + for(int i = 0; i < order; i++) + simplified[i] = simplified[reduced[i]]; + + // fprintf(stderr, "Number of double cosets: %d\n\n", ncosets); + + // simplified_graph = (node_t*) malloc(ncosets*sizeof(node_t)); + seen = (int*) malloc(ncosets*sizeof(int)); + edgelists_used = 0; + + // step 5: set up nodes from minima + current = 0; + for(int i = 0; i < order; i++) + if(reduced[i] == i) { // is minimum + memcpy(simplified_graph[simplified[i]].word, full_graph[i].word, full_graph[i].wordlength*sizeof(int)); + simplified_graph[simplified[i]].wordlength = full_graph[i].wordlength; + simplified_graph[simplified[i]].opposite = simplified[full_graph[i].opposite]; + simplified_graph[simplified[i]].bruhat_lower = (edgelist_t*)0; + simplified_graph[simplified[i]].bruhat_higher = (edgelist_t*)0; + for(int j = 0; j < rank; j++) { + simplified_graph[simplified[i]].left[j] = simplified[full_graph[i].left[j]]; + simplified_graph[simplified[i]].right[j] = simplified[full_graph[i].right[j]]; + } + } + + // step 6: find order relations + for(int i = 0; i < order; i++) { + edge = full_graph[i].bruhat_lower; + while(edge) { + int this = simplified[i]; + int that = simplified[edge->to]; + if(this != that) { + // found something + if(!edgelist_contains(simplified_graph[this].bruhat_lower, that)) + simplified_graph[this].bruhat_lower = edgelist_add(simplified_graph[this].bruhat_lower, that, edgelists_lower, &edgelists_used); + ERROR(simplified_graph[this].wordlength <= simplified_graph[that].wordlength, "The order assumption is being violated!\n"); + } + edge = edge->next; + } + } + + // step 7: remove redundant edges + for(int i = 0; i < ncosets; i++) { + memset(seen, 0, ncosets*sizeof(int)); + queue_init(&queue); + + for(int len = 1; len <= simplified_graph[i].wordlength; len++) { + edge = simplified_graph[i].bruhat_lower; + previous = (edgelist_t*)0; + + while(edge) { + // only look at edges of this length now + if(simplified_graph[i].wordlength - simplified_graph[edge->to].wordlength != len) { + // we only consider edges of length len in this pass + previous = edge; + } else if(seen[edge->to]) { + // this edge is redundant, remove it + // fprintf(stderr, "removing edge from %d to %d\n", i, edge->to); + if(previous) + previous->next = edge->next; + else + simplified_graph[i].bruhat_lower = edge->next; + } else { + // this edge was not redundant, add to seen + previous = edge; + seen[edge->to] = 1; + queue_put(&queue, edge->to); + } + edge = edge->next; + } + + // calculate transitive closure of seen nodes + while((current = queue_get(&queue)) != -1) { + edge = simplified_graph[current].bruhat_lower; + while(edge) { + if(!seen[edge->to]) { + seen[edge->to] = 1; + queue_put(&queue, edge->to); + } + edge = edge->next; + } + } + } + } + + // step 8: revert order + for(int i = 0; i < ncosets; i++) { + edge = simplified_graph[i].bruhat_lower; + while(edge) { + simplified_graph[edge->to].bruhat_higher = + edgelist_add(simplified_graph[edge->to].bruhat_higher, + i, edgelists_higher, &edgelists_used); + edge = edge->next; + } + } + + // output as graphviz dot file + /* + fprintf(stdout, "difull_graph test123 {\n"); + for(int i = 0; i < ncosets; i++) { + edge = simplified_graph[i].bruhat_lower; + while(edge) { + fprintf(stdout, "%s -> %s;\n", + alphabetize(simplified_graph[i].word, simplified_graph[i].wordlength, alphabet, buffer), + alphabetize(simplified_graph[edge->to].word, simplified_graph[edge->to].wordlength, alphabet, buffer2)); + + edge = edge->next; + } + } + fprintf(stdout, "}\n"); */ + + // some output + /* for(int i = 0; i < ncosets; i++) + fprintf(stderr, "%s <=> %s\n", simplified_graph[i].wordlength == 0 ? "1" : alphabetize(simplified_graph[i].word, simplified_graph[i].wordlength, alphabet, buffer), simplified_graph[simplified_graph[i].opposite].wordlength == 0 ? "1" : alphabetize(simplified_graph[simplified_graph[i].opposite].word, simplified_graph[simplified_graph[i].opposite].wordlength, alphabet, buffer2)); */ + + // fprintf(stderr, "\nAdded %d edges.\n\n", edgelists_used); + + free(seen); + free(reduced); + free(group); + free(simplified); + graph_free(type, full_graph); + + return ncosets; +} + +node_t *graph_alloc(semisimple_type_t type) +{ + int rank = coxeter_rank(type); + int order = coxeter_order(type); + int hyperplanes = coxeter_hyperplanes(type); + + node_t *graph = (node_t*)malloc(order*sizeof(node_t)); + int *left = (int*)malloc(order*rank*sizeof(int)); + int *right = (int*)malloc(order*rank*sizeof(int)); + edgelist_t *edgelists = (edgelist_t*)malloc(order*hyperplanes*sizeof(edgelist_t)); + int *words = (int*)malloc(order*hyperplanes*sizeof(int)); + + for(int i = 0; i < order; i++) { + graph[i].left = &left[rank*i]; + graph[i].right = &right[rank*i]; + graph[i].word = &words[hyperplanes*i]; + } + + graph[0].bruhat_higher = edgelists; + + return graph; +} + +void graph_free(semisimple_type_t type, node_t *graph) +{ + free(graph[0].left); + free(graph[0].right); + free(graph[0].word); + + int order = coxeter_order(type); + + // find the head of all edgelists by just taking the one having the lowest address + edgelist_t *edgelists = graph[0].bruhat_lower; + for(int i = 0; i < order; i++) { + if(graph[i].bruhat_lower < edgelists && graph[i].bruhat_lower != 0) + edgelists = graph[i].bruhat_lower; + if(graph[i].bruhat_higher < edgelists && graph[i].bruhat_higher != 0) + edgelists = graph[i].bruhat_higher; + } + free(edgelists); +} + /*********************************** THE ACTUAL ENUMERATION ****************************************/ typedef struct { @@ -401,7 +651,6 @@ static long enumerate_tree(const enumeration_info_t *info, signed char *level, i long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile) { - // int rank, order; signed char *level; long count = 0; enumeration_info_t info; diff --git a/thickenings.h b/thickenings.h index f2f813a..7c6b5d7 100644 --- a/thickenings.h +++ b/thickenings.h @@ -25,10 +25,22 @@ typedef struct { int is_hyperplane_reflection; // boolean value } node_t; +// printing functions char *alphabetize(int *word, int len, const char *alphabet, char *buffer); void print_thickening(int rank, int order, const signed char *thickening, int level, const char *alphabet, FILE *f); -static int compare_wordlength(const void *a, const void *b, void *gr); -void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists_pointer, int **words_pointer); + +// generating the graph of the bruhat order +void prepare_graph(semisimple_type_t type, node_t *graph); +int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigned long right, node_t *simplified_graph); + +// enumerate balanced thickenings long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile); +node_t *graph_alloc(semisimple_type_t type); +void graph_free(semisimple_type_t type, node_t *graph); + +// various helper functions +static int compare_wordlength(const void *a, const void *b, void *gr); +static int edgelist_contains(edgelist_t *list, int x); +static edgelist_t *edgelist_add(edgelist_t *list, int new, edgelist_t *storage, int *storage_index); #endif