diff --git a/Makefile b/Makefile index 467337e..3620e52 100644 --- a/Makefile +++ b/Makefile @@ -1,8 +1,8 @@ HEADERS=weyl.h thickenings.h queue.h bitvec.h -SPECIAL_OPTIONS=-O0 -g +#SPECIAL_OPTIONS=-O0 -g -D_DEBUG #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline -#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline +SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline OPTIONS=-m64 -march=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS) diff --git a/bitvec.h b/bitvec.h index 0494069..c888d32 100644 --- a/bitvec.h +++ b/bitvec.h @@ -107,6 +107,14 @@ static inline void bv_intersection(const bitvec_t *x, const bitvec_t *y, bitvec_ } } +static inline void bv_difference(const bitvec_t *x, const bitvec_t *y, bitvec_t *result) +{ + int i; + for (i=0; i < BV_QWORD_RANK; i++) { + result->v[i] = x->v[i] & ~y->v[i]; + } +} + static inline int bv_disjoint(const bitvec_t *x, const bitvec_t *y) { for(int i = 0; i < BV_QWORD_RANK; i++) diff --git a/generate.c b/generate.c index 3786ec5..58cbc9f 100644 --- a/generate.c +++ b/generate.c @@ -5,15 +5,132 @@ #include #include +char stringbuffer[100]; +char stringbuffer2[100]; + +typedef struct { + node_t *graph; + int cosets; + int rank; + int order; + int hyperplanes; + semisimple_type_t type; + unsigned long left_invariance; + unsigned long right_invariance; + const char *alphabet; + int *buffer; +} info_t; + +int shorten(int i, unsigned long left, unsigned long right, node_t *graph, int rank) +{ + int other, shorter = i; + do { + i = shorter; + for(int j = 0; j < rank; j++) { + other = graph[shorter].left[j]; + if(left & (1 << j) && + graph[other].wordlength < graph[shorter].wordlength) + shorter = other; + other = graph[shorter].right[j]; + if(right & (1 << j) && + graph[other].wordlength < graph[shorter].wordlength) + shorter = other; + } + } while(shorter != i); + + return shorter; +} + void balanced_thickening_callback(const bitvec_t *pos, int size, void *data) { static long totcount = 0; + if(data) { + info_t *info = (info_t*)data; + + unsigned long right_invariance = FIRSTBITS(info->rank); + unsigned long left_invariance = FIRSTBITS(info->rank); + + int bit1, bit2left, bit2right, left, right; + + for(int i = 0; i < size; i++) { + bit1 = i < size/2 ? bv_get_bit(pos, i) : !bv_get_bit(pos, size - 1 - i); + for(int j = 0; j < info->rank; j++) { + left = info->graph[i].left[j]; + right = info->graph[i].right[j]; + bit2left = left < size/2 ? bv_get_bit(pos, left) : !bv_get_bit(pos, size - 1 - left); + bit2right = right < size/2 ? bv_get_bit(pos, right) : !bv_get_bit(pos, size - 1 - right); + if(bit1 != bit2left) + left_invariance &= ~BIT(j); + if(bit1 != bit2right) + right_invariance &= ~BIT(j); + } + } + + printf("left: "); + for(int j = 0; j < info->rank; j++) + printf("%c", left_invariance & (1 << j) ? info->alphabet[j] : ' '); + printf(" right: "); + for(int j = 0; j < info->rank; j++) + printf("%c", right_invariance & (1 << j) ? info->alphabet[j] : ' '); + + if(info->buffer) { + printf(" generators:"); + queue_t queue; + int current, left, right, shortest; + int *buffer = info->buffer; + + for(int i = 0; i < size/2; i++) { + buffer[i] = bv_get_bit(pos, i); + buffer[size-1-i] = !buffer[i]; + } + for(int i = size-1; i >= 0; i--) { + if(buffer[i]) { + int shortest = shorten(i, left_invariance, right_invariance, info->graph, info-> rank); + printf(" %s", alphabetize(info->graph[shortest].word, + info->graph[shortest].wordlength, + info->alphabet, + stringbuffer)); + buffer[i] = 0; + queue_init(&queue); + queue_put(&queue, i); + while((current = queue_get(&queue)) != -1) { + for(edgelist_t *edge = info->graph[current].bruhat_lower; edge != (edgelist_t*)0; edge = edge->next) { + if(buffer[edge->to]) { + buffer[edge->to] = 0; + queue_put(&queue, edge->to); + } + } + for(int j = 0; j < info->rank; j++) { + left = info->graph[current].left[j]; + if(left_invariance & (1 << j) && + info->graph[left].wordlength < info->graph[current].wordlength && + buffer[left]) { + buffer[left] = 0; + queue_put(&queue, left); + } + right = info->graph[current].left[j]; + if(right_invariance & (1 << j) && + info->graph[right].wordlength < info->graph[current].wordlength && + buffer[right]) { + buffer[right] = 0; + queue_put(&queue, right); + } + } + } + } + } + } + + printf("\n"); + } + + /* if((++totcount) % 100000000 == 0) { fprintf(stderr, "Found balanced ideal: "); bv_print(stderr, pos, size/2); fprintf(stderr, "\n"); - } + } */ } int main(int argc, const char *argv[]) @@ -57,8 +174,6 @@ int main(int argc, const char *argv[]) right_invariance |= (1 << (argv[type.n + 2][i] - 'a')); } - ERROR(strlen(alphabet) < rank, "The alphabet has too few letters\n"); - // generate graph graph = graph_alloc(type); @@ -68,48 +183,123 @@ int main(int argc, const char *argv[]) // print stuff + int output_level = 2; + if(getenv("OUTPUT_LEVEL")) + output_level = atoi(getenv("OUTPUT_LEVEL")); + rank = weyl_rank(type); // number of simple roots order = weyl_order(type); // number of Weyl group elements hyperplanes = weyl_hyperplanes(type); // number of positive roots - fprintf(stderr, "Rank: %d\tOrder: %d\tPositive Roots: %d\tCosets: %d\n", rank, order, hyperplanes, cosets); - fprintf(stderr, "\n"); + if(output_level >= 1) { + printf("Poset: "); + if(left_invariance) { + printf("<"); + for(int j = 0; j < rank; j++) + if(left_invariance & BIT(j)) + fputc(alphabet[j], stdout); + printf("> \\ "); + } - /* - 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) { - fprintf(stderr, "\n%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); - wl = graph[i].wordlength; - } else - fprintf(stderr, "%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + for(int i = 0; i < type.n; i++) + printf("%s%c%d", i == 0 ? "" : " x ", type.factors[i].series, type.factors[i].rank); + + if(right_invariance) { + printf(" / <"); + for(int j = 0; j < rank; j++) + if(right_invariance & BIT(j)) + fputc(alphabet[j], stdout); + printf(">"); + } + fprintf(stdout, "\n"); + + fprintf(stdout, "Rank: %d\tOrder: %d\tPositive Roots: %d\tCosets: %d\n\n", rank, order, hyperplanes, cosets); + } + + if(output_level >= 3) { + fprintf(stdout, "Shortest coset representatives: \n"); + for(int i = 0, wl = 0; i < cosets; i++) { + if(i == 0) { + fprintf(stdout, "1"); + } else if(graph[i].wordlength > wl) { + fprintf(stdout, "\n%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + wl = graph[i].wordlength; + } else + fprintf(stdout, "%s ", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + } + fprintf(stdout, "\n\n"); + } + + if(output_level >= 4) { + edgelist_t *edge; + fprintf(stdout, "Bruhat order in graphviz format:\n"); + fprintf(stdout, "digraph test123 {\n"); + for(int i = 0; i < cosets; i++) { + edge = graph[i].bruhat_lower; + while(edge) { + fprintf(stdout, "%s -> %s;\n", + alphabetize(graph[i].word, graph[i].wordlength, alphabet, stringbuffer), + alphabetize(graph[edge->to].word, graph[edge->to].wordlength, alphabet, stringbuffer2)); + + edge = edge->next; + } + } + fprintf(stdout, "}\n\n"); + } + + if(output_level >= 4) { + fprintf(stdout, "Opposites:\n"); + for(int i = 0; i < cosets; i++) { + fprintf(stdout, "%s <-> %s\n", + alphabetize(graph[i].word, graph[i].wordlength, alphabet, stringbuffer), + alphabetize(graph[graph[i].opposite].word, graph[graph[i].opposite].wordlength, alphabet, stringbuffer2)); + + } + fprintf(stdout, "\n"); } - fprintf(stderr, "\n\n"); - */ fixpoints = 0; for(int i = 0; i < cosets; i++) if(graph[i].opposite == i) { - if(fixpoints == 0) - fprintf(stderr, "No thickenings since the longest element fixes the following cosets: %s", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); - else - fprintf(stderr, " %s", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + if(output_level >= 1) { + if(fixpoints == 0) + fprintf(stdout, "No thickenings since the longest element fixes the following cosets: %s", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + else + fprintf(stdout, " %s", alphabetize(graph[i].word, graph[i].wordlength, alphabet, string_buffer1)); + } fixpoints++; } + if(output_level >= 1 && fixpoints) + fprintf(stdout, "\n\n"); - if(fixpoints > 0) { - fprintf(stderr, "\n\n"); - } else { - fwrite(&type.n, sizeof(int), 1, stdout); - fwrite(type.factors, sizeof(simple_type_t), type.n, stdout); - fwrite(&left_invariance, sizeof(unsigned long), type.n, stdout); - fwrite(&right_invariance, sizeof(unsigned long), type.n, stdout); - long count = enumerate_balanced_thickenings(graph, cosets, balanced_thickening_callback, (void*)0); + if(!fixpoints) { + int *buffer = (int*)malloc(cosets*sizeof(int)); - fprintf(stderr, "Found %ld balanced thickenings\n\n", count); + info_t info; + info.graph = graph; + info.cosets = cosets; + info.rank = rank; + info.order = order; + info.hyperplanes = hyperplanes; + info.type = type; + info.left_invariance = left_invariance; + info.right_invariance = right_invariance; + info.alphabet = alphabet; + info.buffer = buffer; + + long count; + if(output_level >= 2) { + fprintf(stdout, "Balanced ideals:\n", count); + count = enumerate_balanced_thickenings(graph, cosets, balanced_thickening_callback, &info); + fprintf(stdout, "\n", count); + } + else + count = enumerate_balanced_thickenings(graph, cosets, 0, 0); + + + if(output_level >= 1) + fprintf(stdout, "Found %ld balanced ideal%s\n", count, count == 1 ? "" : "s"); } graph_free(type, graph); diff --git a/queue.h b/queue.h index a19cc31..ca9b9cc 100644 --- a/queue.h +++ b/queue.h @@ -7,6 +7,11 @@ #define QUEUE_SIZE 5000 #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; diff --git a/thickenings.c b/thickenings.c index 6f60018..9208186 100644 --- a/thickenings.c +++ b/thickenings.c @@ -91,12 +91,10 @@ void prepare_graph(semisimple_type_t type, node_t *graph) graph[i].is_hyperplane_reflection = 0; } - // get coxeter graph + LOG("Generate Weyl group.\n"); weyl_generate(type, graph_data); - fprintf(stderr, "Weyl group generated.\n"); - for(int i = 0; i < order; i++) for(int j = 0; j < rank; j++) { graph_unsorted[i].left = graph_data[i].left; @@ -105,6 +103,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph) // find wordlengths + LOG("Determine word lengths.\n"); + graph_unsorted[0].wordlength = 0; queue_init(&queue); queue_put(&queue, 0); @@ -118,9 +118,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Wordlengths calculated.\n"); - - // sort by wordlength + LOG("Sort by wordlength.\n"); for(int i = 0; i < order; i++) ordering[i] = i; @@ -135,9 +133,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) graph[i].left[j] = reverse_ordering[graph_unsorted[ordering[i]].left[j]]; // rewrite references } - fprintf(stderr, "Sorted by wordlength.\n"); - - // find words + LOG("Find shortest words.\n"); for(int i = 0; i < order; i++) memset(graph[i].word, 0, hyperplanes*sizeof(int)); @@ -154,9 +150,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Shortest words found.\n"); - - // generate right edges + LOG("Generate right edges.\n"); for(int i = 0; i < order; i++) { for(int j = 0; j < rank; j++) { @@ -168,9 +162,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Right edges generated.\n"); - - // find opposites + LOG("Find opposites.\n"); node_t *longest = &graph[order-1]; for(int i = 0; i < order; i++) { @@ -180,9 +172,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) graph[i].opposite = current; } - fprintf(stderr, "Opposites found.\n"); - - // enumerate hyperplanes + LOG("Enumerate hyperplanes.\n"); hyperplane_count = 0; for(int i = 0; i < order; i++) { @@ -203,9 +193,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Hyperplanes enumerated.\n"); - - // generate folding order + LOG("Determine Bruhat order.\n"); edgelist_count = 0; for(int i = 0; i < order; i++) { @@ -229,9 +217,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Bruhat order generated.\n"); - - // remove redundant edges + LOG("Perform transitive reduction.\n"); for(int i = 0; i < order; i++) { memset(seen, 0, order*sizeof(int)); @@ -272,9 +258,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Redundant edges removed.\n"); - - // reverse folding order + LOG("Revert Bruhat order.\n"); edgelist_count = 0; for(int i = 0; i < order; i++) { @@ -288,7 +272,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph) } } - fprintf(stderr, "Bruhat order reversed.\n"); + LOG("Sort opposites.\n"); // additional sorting step to force opposite property (opposite of j is at n - j - 1) @@ -373,7 +357,7 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne full_graph = graph_alloc(type); prepare_graph(type, full_graph); - fprintf(stderr, "Full graph generated.\n"); + LOG("Full graph generated.\n"); // initialize stuff @@ -385,6 +369,8 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne reduced[i] = i; } + LOG("Group by double coset.\n"); + // step 1: group for(int i = 0; i < order; i++) { if(group[i] != -1) @@ -406,6 +392,8 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne } } + LOG("Find minimal length elements.\n"); + // step 2: find minimum for(int i = 0; i < order; i++) if(full_graph[i].wordlength < full_graph[reduced[group[i]]].wordlength) @@ -424,12 +412,11 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne 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; + LOG("Copy minimal elements.\n"); + // step 5: set up nodes from minima current = 0; for(int i = 0; i < order; i++) @@ -445,6 +432,8 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne } } + LOG("Find induced order.\n"); + // step 6: find order relations for(int i = 0; i < order; i++) { edge = full_graph[i].bruhat_lower; @@ -461,6 +450,8 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne } } + LOG("Perform transitive reduction.\n"); + // step 7: remove redundant edges for(int i = 0; i < ncosets; i++) { memset(seen, 0, ncosets*sizeof(int)); @@ -477,7 +468,6 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne 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 @@ -505,6 +495,8 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne } } + LOG("Revert order.\n"); + // step 8: revert order edgelists_used = 0; for(int i = 0; i < ncosets; i++) { @@ -517,33 +509,57 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne } } - // 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)); + LOG("Sort opposites.\n"); - edge = edge->next; + int *ordering = (int*)malloc(ncosets*sizeof(int)); + int *reverse_ordering = (int*)malloc(ncosets*sizeof(int)); + node_t *unsorted = (node_t*)malloc(ncosets*sizeof(node_t)); + int opp, pos; + + pos = 0; + for(int i = 0; i < ncosets; i++) { // first all the pairs + opp = simplified_graph[i].opposite; + if(opp > i) { // first occurrence of this pair + ordering[pos] = i; + ordering[ncosets-1-pos] = opp; + reverse_ordering[i] = pos; + reverse_ordering[opp] = ncosets-1-pos; + pos++; } } - fprintf(stdout, "}\n"); */ + for(int i = 0; i < ncosets; i++) // and finally the self-opposites + if(simplified_graph[i].opposite == i) { + ordering[pos] = i; + reverse_ordering[i] = pos; + pos++; + } - // 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); + // now really do it + memcpy(unsorted, simplified_graph, ncosets*sizeof(node_t)); + for(int i = 0; i < ncosets; i++) { + simplified_graph[i] = unsorted[ordering[i]]; + simplified_graph[i].opposite = reverse_ordering[simplified_graph[i].opposite]; + for(edgelist_t *edge = simplified_graph[i].bruhat_lower; edge != (edgelist_t*)0; edge = edge->next) + edge->to = reverse_ordering[edge->to]; + for(edgelist_t *edge = simplified_graph[i].bruhat_higher; edge != (edgelist_t*)0; edge = edge->next) + edge->to = reverse_ordering[edge->to]; + for(int j = 0; j < rank; j++) { + simplified_graph[i].left[j] = reverse_ordering[simplified_graph[i].left[j]]; + simplified_graph[i].right[j] = reverse_ordering[simplified_graph[i].right[j]]; + } + } + free(ordering); + free(reverse_ordering); + free(unsorted); free(seen); free(reduced); free(group); free(simplified); graph_free(type, full_graph); + LOG("Simplified graph generated.\n"); + return ncosets; } @@ -637,6 +653,11 @@ static long enumerate_tree(const enumeration_info_t *info, const bitvec_t *pos, // everything before next_neg which was unknown should be set to positive; to speed this up, we can start with already_known bv_set_range_except(&newpos, neg, already_known, next_neg); + #ifdef _DEBUG + bv_print_nice(stderr, &newpos, &newneg, -1, info->size/2); + fprintf(stderr, "\n"); + #endif + // check if this leads to any conflicts (equivalently, violates slimness) if(!bv_disjoint(&newpos, &newneg)) return 0; @@ -728,12 +749,15 @@ long enumerate_balanced_thickenings(node_t *graph, int size, void (*callback) (c } free(principal); - /* + // output principal ideals +#ifdef _DEBUG for(int i = 0; i < info.size; i++) { - fprintf(stderr, "%d: ", i); + fprintf(stderr, "%2d: ", i); bv_print_nice(stderr, &info.principal_pos[i], &info.principal_neg[i], -1, info.size/2); fprintf(stderr, "\n"); - } */ + } + fprintf(stderr,"\n"); +#endif // enumerate balanced ideals bitvec_t pos, neg; diff --git a/weyl.c b/weyl.c index 2e05ec3..b595829 100644 --- a/weyl.c +++ b/weyl.c @@ -283,10 +283,10 @@ void weyl_cartan_matrix(semisimple_type_t type, int *m) } break; - case 'B': // not sure at all about the order of B and C + case 'B': // not sure at all about the order of B and C if(type.factors[k].rank >= 2) { - A[0][1] = -1; - A[1][0] = -2; + A[0][1] = -2; + A[1][0] = -1; } for(int i = 2; i < type.factors[k].rank; i++) { A[i][i-1] = -1; @@ -296,8 +296,8 @@ void weyl_cartan_matrix(semisimple_type_t type, int *m) case 'C': if(type.factors[k].rank >= 2) { - A[0][1] = -2; - A[1][0] = -1; + A[0][1] = -1; + A[1][0] = -2; } for(int i = 2; i < type.factors[k].rank; i++) { A[i][i-1] = -1;