diff --git a/Makefile b/Makefile index f0ed9e0..7b6c37c 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,5 @@ HEADERS=coxeter.h thickenings.h queue.h -OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99 -pg +OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99 all: generate process diff --git a/generate.c b/generate.c index cf90399..6e4c87f 100644 --- a/generate.c +++ b/generate.c @@ -72,7 +72,7 @@ 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, alphabet, stdout); + long count = enumerate_balanced_thickenings(type, graph, order, alphabet, stdout); fprintf(stderr, "\n"); fprintf(stderr, "Found %ld balanced thickenings\n\n", count); diff --git a/test.c b/test.c new file mode 100644 index 0000000..b9600ee --- /dev/null +++ b/test.c @@ -0,0 +1,257 @@ +#include +#include + +#include "thickenings.h" +#include "queue.h" + +#define SWAP(t, a, b) do {t tmp = a; a = b; b = tmp;} while(0) + +int edgelist_contains(edgelist_t *list, int needle) { + while(list) { + if(list->to == needle) + return 1; + list = list->next; + } + return 0; +} + +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 main(int argc, const char *argv[]) +{ + semisimple_type_t type; + node_t *graph; + int *leftbuf, *rightbuf; + edgelist_t *edgelists; + edgelist_t *edgelists_simplified; + int edgelists_simplified_used; + int *words; + int rank, order, max_wordlength; + int *reduced, *group, *simplified; + int *seen; + int current; + edgelist_t *edge, *previous; + queue_t queue; + + char alphabet[] = "abcdefghijklmnopqrstuvwxyz"; + char buffer[1024], buffer2[1024]; + + int ncosets; + node_t *simplified_graph; + + // left and right invariances as bitmasks + // int left = ~(1 << (atoi(argv[1]) - atoi(argv[3]))); + // int right = ~(1 << (atoi(argv[1]) - atoi(argv[2]))); + int left = atoi(argv[2]); + int right = atoi(argv[3]); + + 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]); + + rank = coxeter_rank(type); + order = coxeter_order(type); + graph = (node_t*)malloc(order*sizeof(node_t)); + leftbuf = (int*)malloc(rank*order*sizeof(int)); + rightbuf = (int*)malloc(rank*order*sizeof(int)); + for(int i = 0; i < order; i++) { + graph[i].left = &leftbuf[i*rank]; + graph[i].right = &rightbuf[i*rank]; + } + + prepare_graph(type, graph, &edgelists, &words); + + 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, graph[current].left[j]); + if(right & (1 << j)) + queue_put(&queue, graph[current].right[j]); + } + } + } + + // step 2: find minimum + for(int i = 0; i < order; i++) + if(graph[i].wordlength < 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]]; + + // count 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); + + max_wordlength = coxeter_hyperplanes(type); + simplified_graph = (node_t*) malloc(ncosets*sizeof(node_t)); + edgelists_simplified = (edgelist_t*) malloc(2*max_wordlength*order*sizeof(edgelist_t)); + seen = (int*) malloc(ncosets*sizeof(int)); + edgelists_simplified_used = 0; + + // copy minima + current = 0; + for(int i = 0; i < order; i++) + if(reduced[i] == i) { // is minimum + simplified_graph[simplified[i]].word = graph[i].word; + simplified_graph[simplified[i]].wordlength = graph[i].wordlength; + simplified_graph[simplified[i]].opposite = simplified[graph[i].opposite]; + simplified_graph[simplified[i]].bruhat_lower = (edgelist_t*)0; + simplified_graph[simplified[i]].bruhat_higher = (edgelist_t*)0; + } + + // 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)); + + // find order relations + for(int i = 0; i < order; i++) { + edge = 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_simplified, &edgelists_simplified_used); + // if(!edgelist_contains(simplified_graph[that].bruhat_higher, this)) + // simplified_graph[that].bruhat_higher = edgelist_add(simplified_graph[that].bruhat_higher, this, edgelists_simplified, &edgelists_simplified_used); + ERROR(simplified_graph[this].wordlength <= simplified_graph[that].wordlength, "The order assumption is being violated!\n"); + } + edge = edge->next; + } + } + + fprintf(stderr, "\nAdded %d edges.\n\n", edgelists_simplified_used); + + // 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; + } + } + } + } + + // reverse 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_simplified, &edgelists_simplified_used); + edge = edge->next; + } + } + + + fprintf(stdout, "digraph 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"); + + long nthickenings = enumerate_balanced_thickenings(type, simplified_graph, ncosets, alphabet, stdout); + + fprintf(stderr, "Found %ld balanced thickenings.\n", nthickenings); + + /* for(int i = 0; i < order; i++) + printf("%s <= %s\n", + reduced[i] == 0 ? "1" : alphabetize(graph[reduced[i]].word, graph[reduced[i]].wordlength, alphabet, buffer2), + i == 0 ? "1" : alphabetize(graph[i].word, graph[i].wordlength, alphabet, buffer)); */ + + + free(seen); + free(simplified_graph); + free(edgelists_simplified); + free(type.factors); + free(graph); + free(edgelists); + free(words); + free(reduced); + free(group); + free(simplified); + free(leftbuf); + free(rightbuf); + + return 0; +} diff --git a/thickenings.c b/thickenings.c index 1ef6f19..f468733 100644 --- a/thickenings.c +++ b/thickenings.c @@ -289,6 +289,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph, edgelist_t **edgelists typedef struct { int rank; int order; + int size; // the size of the graph; this can vary from the order if we take quotients beforehand const node_t *graph; int printstep; const char *alphabet; @@ -307,7 +308,7 @@ static int transitive_closure(const enumeration_info_t *info, signed char *level level[info->graph[head].opposite] = -current_level; queue_put(&queue, head); - for(int i = head + 1; level[i] != HEAD_MARKER && i < info->order; i++) { // everything which is right to the head and empty will not get marked in this or higher levels, so we can mark its opposite + for(int i = head + 1; level[i] != HEAD_MARKER && i < info->size; i++) { // everything which is right to the head and empty will not get marked in this or higher levels, so we can mark its opposite if(level[i] == current_level) { is_slim = 0; break; @@ -342,11 +343,11 @@ static inline void output_thickening(const enumeration_info_t *info, signed char { // if printstep is set accordingly, write state to stderr if(is_slim && is_fat && info->printstep > 0 && (count + 1) % info->printstep == 0) { - print_thickening(info->rank, info->order, level, current_level, info->alphabet, stderr); + print_thickening(info->rank, info->size, level, current_level, info->alphabet, stderr); fprintf(stderr, "\n"); } else if(info->printstep < 0) { - print_thickening(info->rank, info->order, level, current_level - !is_slim, info->alphabet, stderr); + print_thickening(info->rank, info->size, level, current_level - !is_slim, info->alphabet, stderr); fprintf(stderr, " "); if(is_slim) { fprintf(stderr, "S"); @@ -378,7 +379,7 @@ static long enumerate_tree(const enumeration_info_t *info, signed char *level, i if(is_fat) { count++; - fwrite(level, sizeof(signed char), info->order, info->outfile); + fwrite(level, sizeof(signed char), info->size, info->outfile); } else { for(int i = head - 1; i >= 0; i--) if(level[i] == 0) @@ -388,14 +389,14 @@ static long enumerate_tree(const enumeration_info_t *info, signed char *level, i // clean up level[head] = 0; - for(int i = 0; i < info->order; i++) + for(int i = 0; i < info->size; i++) if(level[i] >= current_level && level[i] != HEAD_MARKER || level[i] <= -current_level) level[i] = 0; return count; } -long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const char *alphabet, FILE *outfile) +long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile) { // int rank, order; signed char *level; @@ -406,6 +407,7 @@ long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const info.rank = coxeter_rank(type); info.order = coxeter_order(type); + info.size = size; info.graph = graph; info.alphabet = (char*)alphabet; info.outfile = outfile; @@ -414,10 +416,16 @@ long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const if(getenv("PRINTSTEP")) info.printstep = atoi(getenv("PRINTSTEP")); - level = (signed char*)malloc(info.order*sizeof(int)); - memset(level, 0, info.order*sizeof(int)); + // the algorithm only works if the opposition pairing does not stabilize any element + // if this happens, there can be no balanced thickenings + for(int i = 0; i < info.size; i++) + if(graph[i].opposite == i) + return 0; - for(int i = info.order - 1; i >= 0; i--) + level = (signed char*)malloc(info.size*sizeof(int)); + memset(level, 0, info.size*sizeof(int)); + + for(int i = info.size - 1; i >= 0; i--) count += enumerate_tree(&info, level, 1, i); free(level); diff --git a/thickenings.h b/thickenings.h index e762325..f2f813a 100644 --- a/thickenings.h +++ b/thickenings.h @@ -29,6 +29,6 @@ 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); -long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, const char *alphabet, FILE *outfile); +long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int size, const char *alphabet, FILE *outfile); #endif