From efd8e621ea39cde83ab5005ca37099c061c04d65 Mon Sep 17 00:00:00 2001
From: Florian Stecker <mail@florianstecker.de>
Date: Fri, 27 Jan 2017 20:48:44 +0100
Subject: [PATCH] Major rewrite

---
 Makefile      |  12 +-
 generate.c    | 143 ++++++------
 thickenings.c | 629 ++------------------------------------------------
 thickenings.h |  39 +---
 weyl.c        | 617 ++++++++++++++++++++++++++++++++++++++-----------
 weyl.h        |  94 ++++++--
 6 files changed, 641 insertions(+), 893 deletions(-)

diff --git a/Makefile b/Makefile
index 3620e52..ad7125d 100644
--- a/Makefile
+++ b/Makefile
@@ -6,20 +6,14 @@ SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 
 OPTIONS=-m64 -march=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
 
-all: generate process
+all: enumerate
 
-generate: generate.o weyl.o thickenings.o
+enumerate: generate.o weyl.o thickenings.o
 	gcc $(OPTIONS) -o generate generate.o thickenings.o weyl.o
 
-process: process.o weyl.o thickenings.o
-	gcc $(OPTIONS) -o process process.o thickenings.o weyl.o
-
 generate.o: generate.c $(HEADERS)
 	gcc $(OPTIONS) -c generate.c
 
-process.o: process.c $(HEADERS)
-	gcc $(OPTIONS) -c process.c
-
 thickenings.o: thickenings.c $(HEADERS)
 	gcc $(OPTIONS) -c thickenings.c
 
@@ -27,4 +21,4 @@ weyl.o: weyl.c $(HEADERS)
 	gcc $(OPTIONS) -c weyl.c
 
 clean:
-	rm -f generate process thickenings.o weyl.o generate.o process.o
+	rm -f generate thickenings.o weyl.o generate.o
diff --git a/generate.c b/generate.c
index c069cd7..cdbdd70 100644
--- a/generate.c
+++ b/generate.c
@@ -9,20 +9,27 @@ char stringbuffer[100];
 char stringbuffer2[100];
 
 typedef struct {
-  node_t *graph;
-  int cosets;
+  doublequotient_t *dq;
   int rank;
   int order;
-  int hyperplanes;
-  semisimple_type_t type;
-  unsigned long left_invariance;
-  unsigned long right_invariance;
-  const char *alphabet;
+  int positive;
   int *buffer;
   int level;
 } info_t;
 
-int shorten(int i, unsigned long left, unsigned long right, node_t *graph, int rank)
+static char* alphabetize(weylgroup_element_t *e, char *str)
+{
+  if(e->wordlength == 0)
+    sprintf(str, "1");
+  else
+    for(int j = 0; j < e->wordlength; j++)
+      sprintf(str, "%c", e->word[j] + 'a');
+
+  return str;
+}
+
+/*
+int shorten(int i, unsigned long left, unsigned long right, doublequotient_t *dq, int rank)
 {
   int other, shorter = i;
   do {
@@ -41,6 +48,7 @@ int shorten(int i, unsigned long left, unsigned long right, node_t *graph, int r
 
   return shorter;
 }
+*/
 
 void balanced_thickening_callback(const bitvec_t *pos, int size, void *data)
 {
@@ -57,8 +65,8 @@ void balanced_thickening_callback(const bitvec_t *pos, int size, void *data)
     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];
+	left = info->dq->cosets[i].min->left[j]->coset - info->dq->cosets;
+	right = info->dq->cosets[i].min->right[j]->coset - info->dq->cosets;
 	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)
@@ -70,15 +78,16 @@ void balanced_thickening_callback(const bitvec_t *pos, int size, void *data)
 
     printf("%4d left: ", totcount++);
     for(int j = 0; j < info->rank; j++)
-      printf("%c", left_invariance & (1 << j) ? info->alphabet[j] : ' ');
+      printf("%c", left_invariance & (1 << j) ? j + 'a' : ' ');
     printf(" right: ");
     for(int j = 0; j < info->rank; j++)
-      printf("%c", right_invariance & (1 << j) ? info->alphabet[j] : ' ');
+      printf("%c", right_invariance & (1 << j) ? j + 'a' : ' ');
 
+    /*
     if(info->buffer) {
       printf(" generators:");
       queue_t queue;
-      int current, left, right, shortest;
+      int cur, left, right;
       int *buffer = info->buffer;
 
       for(int i = 0; i < size/2; i++) {
@@ -87,32 +96,30 @@ void balanced_thickening_callback(const bitvec_t *pos, int size, void *data)
       }
       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));
+	  weylgroup_element_t *shortest = shorten(i, left_invariance, right_invariance, info->dq);
+	  printf(" %s", alphabetize(shortest, 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);
+	  while((cur = queue_get(&queue)) != -1) {
+	    for(doublecoset_list_t *current = info->dq->coset[cur].bruhat_lower; current; current = current->next) {
+	      int idx = current->to - info->dq->cosets;
+	      if(buffer[idx]) {
+		buffer[idx] = 0;
+		queue_put(&queue, idx);
 	      }
 	    }
 	    for(int j = 0; j < info->rank; j++) {
-	      left = info->graph[current].left[j];
+	      left = info->dq->coset[cur].min.left[j] - info->dq->cosets;
 	      if(left_invariance & (1 << j) &&
-		 info->graph[left].wordlength < info->graph[current].wordlength &&
+		 info->graph[left].wordlength < info->graph[cur].wordlength &&
 		 buffer[left]) {
 		buffer[left] = 0;
 		queue_put(&queue, left);
 	      }
-	      right = info->graph[current].left[j];
+	      right = info->graph[cur].left[j];
 	      if(right_invariance & (1 << j) &&
-		 info->graph[right].wordlength < info->graph[current].wordlength &&
+		 info->graph[right].wordlength < info->graph[cur].wordlength &&
 		 buffer[right]) {
 		buffer[right] = 0;
 		queue_put(&queue, right);
@@ -132,6 +139,7 @@ void balanced_thickening_callback(const bitvec_t *pos, int size, void *data)
 	printf("]");
       }
     }
+    */
 
     printf("\n");
   }
@@ -151,12 +159,11 @@ int main(int argc, const char *argv[])
 {
   semisimple_type_t type;
   unsigned long right_invariance, left_invariance;
-  int rank, order, hyperplanes, cosets;
+  int rank, order, positive;
   int fixpoints;
 
-  node_t *graph;
+  doublequotient_t *dq;
 
-  char string_buffer1[1000];
   const char *alphabet = "abcdefghijklmnopqrstuvwxyz";
 
   // read arguments
@@ -190,10 +197,7 @@ int main(int argc, const char *argv[])
 
   // generate graph
 
-  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!\n");
+  dq = weyl_generate_bruhat(type, left_invariance, right_invariance);
 
   // print stuff
 
@@ -203,7 +207,7 @@ int main(int argc, const char *argv[])
 
   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
+  positive = weyl_positive(type);        // number of positive roots
 
   if(output_level >= 1) {
     if(left_invariance) {
@@ -226,59 +230,48 @@ int main(int argc, const char *argv[])
     }
     fprintf(stdout, "\n");
 
-    fprintf(stdout, "Rank: %d\tOrder: %d\tPositive Roots: %d\tCosets: %d\n\n", rank, order, hyperplanes, cosets);
+    fprintf(stdout, "Rank: %d\tOrder: %d\tPositive Roots: %d\tCosets: %d\n\n", rank, order, positive, dq->count);
   }
 
   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));
+    for(int i = 0, wl = 0; i < dq->count; i++) {
+      if(dq->cosets[i].min->wordlength > wl) {
+	printf("\n");
+	wl = dq->cosets[i].min->wordlength;
+      }
+      fprintf(stdout, "\n%s ", alphabetize(dq->cosets[i].min, stringbuffer));
     }
     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) {
+    for(int i = 0; i < dq->count; i++)
+      for(doublecoset_list_t *current = dq->cosets[i].bruhat_lower; current; current = current->next)
 	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;
-      }
-    }
+		alphabetize(dq->cosets[i].min, stringbuffer),
+		alphabetize(current->to->min, stringbuffer2));
     fprintf(stdout, "}\n\n");
   }
 
   if(output_level >= 4) {
     fprintf(stdout, "Opposites:\n");
-    for(int i = 0; i < cosets; i++) {
+    for(int i = 0; i < dq->count; 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));
-
-    }
+	      alphabetize(dq->cosets[i].min, stringbuffer),
+	      alphabetize(dq->cosets[i].opposite->min, stringbuffer2));
     fprintf(stdout, "\n");
   }
 
   fixpoints = 0;
-  for(int i = 0; i < cosets; i++)
-    if(graph[i].opposite == i) {
+  for(int i = 0; i < dq->count; i++)
+    if(dq->cosets[i].opposite == &dq->cosets[i]) {
       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));
+	  fprintf(stdout, "No thickenings since the longest element fixes the following cosets:");
+	fprintf(stdout, " %s", alphabetize(dq->cosets[i].min, stringbuffer));
       }
       fixpoints++;
     }
@@ -287,37 +280,31 @@ int main(int argc, const char *argv[])
 
 
   if(!fixpoints) {
-    int *buffer = (int*)malloc(cosets*sizeof(int));
+    int *buffer = (int*)malloc(dq->count*sizeof(int));
 
     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.dq = dq;
+    info.rank = weyl_rank(type);
+    info.order = weyl_order(type);
+    info.positive = weyl_positive(type);
     info.buffer = buffer;
     info.level = output_level;
 
     long count;
     if(output_level >= 2) {
       fprintf(stdout, "Balanced ideals:\n", count);
-      count = enumerate_balanced_thickenings(graph, cosets, balanced_thickening_callback, &info);
+      count = enumerate_balanced_thickenings(dq, balanced_thickening_callback, &info);
       fprintf(stdout, "\n", count);
     } else {
       long outputcount = 0;
-      count = enumerate_balanced_thickenings(graph, cosets, balanced_thickening_simple_callback, &outputcount);
+      count = enumerate_balanced_thickenings(dq, balanced_thickening_simple_callback, &outputcount);
     }
 
-
     if(output_level >= 1)
       fprintf(stdout, "Found %ld balanced ideal%s\n", count, count == 1 ? "" : "s");
   }
 
-  graph_free(type, graph);
+  weyl_destroy_bruhat(dq);
   free(type.factors);
 
   return 0;
diff --git a/thickenings.c b/thickenings.c
index 7794ded..58d1cd2 100644
--- a/thickenings.c
+++ b/thickenings.c
@@ -8,606 +8,6 @@
 #include "weyl.h"
 #include "queue.h"
 
-char *alphabetize(int *word, int len, const char *alphabet, char *buffer)
-{
-  if(len == 0) {
-    buffer[0] = '1';
-    buffer[1] = 0;
-    return buffer;
-  }
-
-  int i = 0;
-  for(i = 0; i < len; i++)
-    buffer[i] = alphabet[word[i]];
-  buffer[i] = 0;
-
-  return buffer;
-}
-
-void print_thickening(int rank, int order, const signed char *thickening, int upto_level, const char *alphabet, FILE *f)
-{
-  for(int i = 0; i < order; i++) {
-    if(thickening[i] == HEAD_MARKER)
-      fprintf(f, "\e[41;37mx\e[0m");
-    else if(thickening[i] < - upto_level || thickening[i] > upto_level)
-      fprintf(f, " ");
-    else if(thickening[i] < 0 && thickening[i] > -10)
-      fprintf(f, "\e[47;30m%d\e[0m", -thickening[i]);
-    else if(thickening[i] <= -10)
-      fprintf(f, "\e[47;30m+\e[0m");
-    else if(thickening[i] > 0 && thickening[i] < 10)
-      fprintf(f, "\e[40;37m%d\e[0m", thickening[i]);
-    else if(thickening[i] >= 10)
-      fprintf(f, "\e[40;37m+\e[0m");
-    else
-      fprintf(f, " ");
-  }
-
-  fprintf(f, "\e[K");
-}
-
-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;
-}
-
-void prepare_graph(semisimple_type_t type, node_t *graph)
-{
-  queue_t queue;
-
-  edgelist_t *edgelists_lower, *edgelists_higher;
-  int rank, order, hyperplanes;
-  edgelist_t *edge, *previous;
-  int edgelist_count, hyperplane_count;
-  int current;
-
-  weylgroup_element_t *graph_data;
-  node_t *graph_unsorted;
-  int *ordering, *reverse_ordering, *seen;
-
-  // initialize
-
-  rank = weyl_rank(type);
-  order = weyl_order(type);
-  hyperplanes = weyl_hyperplanes(type);
-
-  edgelists_higher = graph[0].bruhat_higher;
-  edgelists_lower = &graph[0].bruhat_higher[order*hyperplanes/2];
-
-  graph_data = weyl_alloc(type);
-  graph_unsorted = (node_t*)malloc(order*sizeof(node_t));
-  ordering = (int*)malloc(order*sizeof(int));
-  reverse_ordering = (int*)malloc(order*sizeof(int));
-  seen = (int*)malloc(order*sizeof(int));
-
-  for(int i = 0; i < order; i++) {
-    graph_unsorted[i].wordlength = INT_MAX;
-    graph[i].bruhat_lower = 0;
-    graph[i].bruhat_higher = 0;
-    graph[i].is_hyperplane_reflection = 0;
-  }
-
-  LOG("Generate Weyl group.\n");
-
-  weyl_generate(type, graph_data);
-
-  for(int i = 0; i < order; i++)
-    for(int j = 0; j < rank; j++) {
-      graph_unsorted[i].left = graph_data[i].left;
-      graph_unsorted[i].id = graph_data[i].id;
-    }
-
-  // find wordlengths
-
-  LOG("Determine word lengths.\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);
-      }
-    }
-  }
-
-  LOG("Sort by wordlength.\n");
-
-  for(int i = 0; i < order; i++)
-    ordering[i] = i;
-  qsort_r(ordering, order, sizeof(int), compare_wordlength, graph_unsorted); // so ordering is a map new index -> old index
-  for(int i = 0; i < order; i++)
-    reverse_ordering[ordering[i]] = i; // reverse_ordering is a map old index -> new index
-  for(int i = 0; i < order; i++) {
-    // we have only set left, wordlength and id so far, so just copy these
-    graph[i].wordlength = graph_unsorted[ordering[i]].wordlength;
-    graph[i].id = graph_unsorted[ordering[i]].id;
-    for(int j = 0; j < rank; j++)
-      graph[i].left[j] = reverse_ordering[graph_unsorted[ordering[i]].left[j]]; // rewrite references
-  }
-
-  LOG("Find shortest words.\n");
-
-  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] == 0) {
-	memcpy(&graph[neighbor].word[1], &graph[current].word[0], graph[current].wordlength*sizeof(int));
-	graph[neighbor].word[0] = i;
-	queue_put(&queue, neighbor);
-      }
-    }
-  }
-
-  LOG("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;
-    }
-  }
-
-  LOG("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;
-  }
-
-  LOG("Enumerate hyperplanes.\n");
-
-  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++;
-      }
-    }
-  }
-
-  LOG("Determine Bruhat order.\n");
-
-  edgelist_count = 0;
-  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_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 {
-	  ERROR(1, "Chambers of equal word lengths should not be folded on each other!\n");
-	}
-      }
-    }
-  }
-
-  LOG("Perform transitive reduction.\n");
-
-  for(int i = 0; i < order; i++) {
-    memset(seen, 0, order*sizeof(int));
-    queue_init(&queue);
-
-    for(int len = 1; len <= graph[i].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(graph[i].wordlength - graph[edge->to].wordlength != len) {
-	  previous = edge;
-	} else if(seen[edge->to]) {
-	  if(previous)
-	    previous->next = edge->next;
-	  else
-	    graph[i].bruhat_lower = edge->next;
-	} else {
-	  previous = edge;
-	  seen[edge->to] = 1;
-	  queue_put(&queue, edge->to);
-	}
-	edge = edge->next;
-      }
-
-      // see which nodes we can reach using only edges up to length len, mark them as seen
-      while((current = queue_get(&queue)) != -1) {
-	edge = graph[current].bruhat_lower;
-	while(edge) {
-	  if(!seen[edge->to]) {
-	    seen[edge->to] = 1;
-	    queue_put(&queue, edge->to);
-	  }
-	  edge = edge->next;
-	}
-      }
-    }
-  }
-
-  LOG("Revert Bruhat order.\n");
-
-  edgelist_count = 0;
-  for(int i = 0; i < order; i++) {
-    edge = graph[i].bruhat_lower;
-    while(edge) {
-      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;
-    }
-  }
-
-  LOG("Sort opposites.\n");
-
-  // additional sorting step to force opposite property (opposite of j is at n - j - 1)
-
-  for(int i = 0; i < order; i++)
-    reverse_ordering[i] = -1;
-  for(int i = 0, j = 0; i < order; i++) {  // i = old index, j = new index
-    if(reverse_ordering[i] == -1) {
-      reverse_ordering[i] = j;
-      ordering[j] = i;
-      reverse_ordering[graph[i].opposite] = order - j - 1;
-      ordering[order - j - 1] = graph[i].opposite;
-      j++;
-    }
-  }
-  memcpy(graph_unsorted, graph, order*sizeof(node_t));
-  for(int i = 0; i < order; i++) {
-    graph[i] = graph_unsorted[ordering[i]];
-    graph[i].opposite = reverse_ordering[graph[i].opposite];
-    for(int j = 0; j < rank; j++) {
-      graph[i].left[j] = reverse_ordering[graph[i].left[j]];
-      graph[i].right[j] = reverse_ordering[graph[i].right[j]];
-    }
-    for(edge = graph[i].bruhat_lower; edge; edge = edge->next)
-      edge->to = reverse_ordering[edge->to];
-    for(edge = graph[i].bruhat_higher; edge; edge = edge->next)
-      edge->to = reverse_ordering[edge->to];
-  }
-
-  weyl_free(graph_data);
-  free(graph_unsorted);
-  free(ordering);
-  free(reverse_ordering);
-  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;
-
-  rank = weyl_rank(type);
-  order = weyl_order(type);
-  hyperplanes = weyl_hyperplanes(type);
-
-  for(int i = 0; i < rank; i++) {
-    int oppi = weyl_opposition(type, i);
-    if(left & BIT(i) && !(left & BIT(oppi)) ||
-       left & BIT(oppi) && !(left & BIT(i)))
-      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);
-
-  LOG("Full graph generated.\n");
-
-  // initialize stuff
-
-  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;
-  }
-
-  LOG("Group by double coset.\n");
-
-  // 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]);
-      }
-    }
-  }
-
-  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)
-      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]];
-
-  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++)
-    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]].id = full_graph[i].id;
-      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]];
-      }
-    }
-
-  LOG("Find induced order.\n");
-
-  // 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;
-    }
-  }
-
-  LOG("Perform transitive reduction.\n");
-
-  // 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
-	  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;
-	}
-      }
-    }
-  }
-
-  LOG("Revert order.\n");
-
-  // step 8: revert order
-  edgelists_used = 0;
-  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;
-    }
-  }
-
-  LOG("Sort opposites.\n");
-
-  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++;
-    }
-  }
-  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++;
-    }
-
-  // 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;
-}
-
-node_t *graph_alloc(semisimple_type_t type)
-{
-  int rank = weyl_rank(type);
-  int order = weyl_order(type);
-  int hyperplanes = weyl_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 = weyl_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 {
   int size;        // the size of the weyl group. We store however only the first size/2 elements
   bitvec_t *principal_pos;
@@ -680,18 +80,15 @@ static long enumerate_tree(const enumeration_info_t *info, const bitvec_t *pos,
     next_next_neg = bv_next_zero(&known, next_next_neg + 1);
   } while(next_next_neg <= info->size/2);
 
-  // multiprocessing stuff
-  //  if(level == 3)
-  //    fprintf(stderr, "%d\n", count);
-
   return count;
 }
 
-void generate_principal_ideals(node_t *graph, int size, bitvec_t *pos, bitvec_t *neg, int *is_slim)
+static void generate_principal_ideals(doublequotient_t *dq, bitvec_t *pos, bitvec_t *neg, int *is_slim)
 {
   queue_t queue;
   int current;
-  edgelist_t *edge;
+  doublecoset_list_t *edge;
+  int size = dq->count;
 
   // generate principal ideals
   int *principal = (int*)malloc(size*sizeof(int));
@@ -701,10 +98,10 @@ void generate_principal_ideals(node_t *graph, int size, bitvec_t *pos, bitvec_t
     queue_init(&queue);
     queue_put(&queue, i);
     while((current = queue_get(&queue)) != -1)
-      for(edge = graph[current].bruhat_lower; edge; edge = edge->next)
-	if(!principal[edge->to]) {
-	  principal[edge->to] = 1;
-	  queue_put(&queue, edge->to);
+      for(edge = dq->cosets[current].bruhat_lower; edge; edge = edge->next)
+	if(!principal[edge->to - dq->cosets]) {
+	  principal[edge->to - dq->cosets] = 1;
+	  queue_put(&queue, edge->to - dq->cosets);
 	}
 
     // copy the first half into bitvectors
@@ -726,7 +123,7 @@ void generate_principal_ideals(node_t *graph, int size, bitvec_t *pos, bitvec_t
       fprintf(stderr, " ids: [0");
       for(int j = 1; j < size; j++)
 	if(principal[j])
-	  fprintf(stderr, ", %d", graph[j].id);
+	  fprintf(stderr, ", %d", dq->cosets[j].min->id);
       fprintf(stderr, "]\n");
     }
 #endif
@@ -757,12 +154,12 @@ void generate_principal_ideals(node_t *graph, int size, bitvec_t *pos, bitvec_t
    returns the number of balanced ideals
 */
 
-long enumerate_balanced_thickenings(node_t *graph, int size, void (*callback) (const bitvec_t *, int, void*), void *callback_data)
+long enumerate_balanced_thickenings(doublequotient_t *dq, void (*callback) (const bitvec_t *, int, void*), void *callback_data)
 {
   long count = 0;
   enumeration_info_t info;
 
-  info.size = size;
+  info.size = dq->count;
   info.callback = callback;
   info.callback_data = callback_data;
   info.principal_pos = (bitvec_t*)malloc(info.size*sizeof(bitvec_t));
@@ -771,15 +168,15 @@ long enumerate_balanced_thickenings(node_t *graph, int size, void (*callback) (c
 
   // 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)
+  for(int i = 0; i < dq->count; i++)
+    if(dq->cosets[i].opposite->min->id == dq->cosets[i].min->id)
       return 0;
 
   // we can only handle bitvectors up to 64*BV_QWORD_RANK bits, but we only store half of the weyl group
   if(info.size > 128*BV_QWORD_RANK)
     return -1;
 
-  generate_principal_ideals(graph, size, info.principal_pos, info.principal_neg, info.principal_is_slim);
+  generate_principal_ideals(dq, info.principal_pos, info.principal_neg, info.principal_is_slim);
 
   // enumerate balanced ideals
   bitvec_t pos, neg;
diff --git a/thickenings.h b/thickenings.h
index 5e59806..b8b8e8e 100644
--- a/thickenings.h
+++ b/thickenings.h
@@ -3,48 +3,11 @@
 
 #define BV_QWORD_RANK 10
 #include "bitvec.h"
-
 #include "weyl.h"
 
 #define DEBUG(msg, ...) do{fprintf(stderr, msg, ##__VA_ARGS__); }while(0)
 
-#define MAX_THICKENINGS 0 // 0 means infinite
-#define HEAD_MARKER 127
-
-typedef struct _edgelist {
-  int to;
-  struct _edgelist *next;
-} edgelist_t;
-
-// describes an element of the Weyl group; only "opposite" and "bruhat_lower" are being used for enumerating thickenings; everything else is just needed for initialization or output
-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
-  weylid_t id;
-} 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);
-
-// generating the graph of the bruhat order
-node_t *graph_alloc(semisimple_type_t type);
-void graph_free(semisimple_type_t type, node_t *graph);
-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);
-
 // enumerating balanced thickenings
-long enumerate_balanced_thickenings(node_t *graph, int size, void (*callback) (const bitvec_t *, int, void*), void *callback_data);
-
-// 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);
+long enumerate_balanced_thickenings(doublequotient_t *dq, void (*callback) (const bitvec_t *, int, void*), void *callback_data);
 
 #endif
diff --git a/weyl.c b/weyl.c
index 30319dc..7c4fbff 100644
--- a/weyl.c
+++ b/weyl.c
@@ -12,15 +12,167 @@ typedef struct {
   int position;
 } weylid_lookup_t;
 
+static void generate_left_and_ids(semisimple_type_t type, weylgroup_element_t *group);
 static int search(const void *key, const void *base, size_t nmem, size_t size, int (*compar) (const void *, const void *, void *), void *arg);
 static int compare_root_vectors(int rank, const int *x, const int *y);
 static int compare_root_vectors_qsort(const void *x, const void *y, void *arg);
+static int compare_weylid(const void *x, const void *y);
 static int compare_weylid_lookup(const void *x, const void *y);
 static int lookup_id(weylid_t id, weylid_lookup_t *list, int len);
 static weylid_t multiply_generator(int s, weylid_t w, const int *simple, const int *mapping, int rank, int positive);
 static void reflect_root_vector(const int *cartan, int rank, int i, int *old, int *new);
+static weylgroup_element_t* apply_word(int *word, int len, weylgroup_element_t *current);
+static weylgroup_element_t* apply_word_reverse(int *word, int len, weylgroup_element_t *current);
 
-/***************** simple helper functions **********************************/
+/******** generate_left_and_ids and a pile of helper functions **************/
+
+static void generate_left_and_ids(semisimple_type_t type, weylgroup_element_t *group)
+{
+  int rank = weyl_rank(type);
+  int order = weyl_order(type);
+  int positive = weyl_positive(type);
+
+  queue_t queue;
+  int current;
+  int roots_known, elements, length_elements, nextids_count;
+  int *cartan_matrix;
+  int *root_vectors;
+  int *vector;
+  int *simple_roots;
+  int *root_mapping;
+  weylid_t *ids, *edges, *nextids;
+  weylid_lookup_t *lookup;
+
+  // allocate temporary stuff
+
+  cartan_matrix =      (int*)malloc(rank*rank      *sizeof(int));
+  root_vectors =       (int*)malloc(2*positive*rank*sizeof(int));
+  vector =             (int*)malloc(rank           *sizeof(int));
+  root_mapping =       (int*)malloc(positive*rank  *sizeof(int));
+  simple_roots =       (int*)malloc(rank           *sizeof(int));
+  ids =           (weylid_t*)malloc(order          *sizeof(weylid_t));
+  edges =         (weylid_t*)malloc(rank*order     *sizeof(weylid_t));
+  nextids =       (weylid_t*)malloc(rank*order     *sizeof(weylid_t));
+  lookup = (weylid_lookup_t*)malloc(order          *sizeof(weylid_lookup_t));
+
+  // get all information on the cartan type
+  LOG("Get Cartan matrix.\n");
+
+  weyl_cartan_matrix(type, cartan_matrix);
+
+  // enumerate roots, first the simple ones, then all others by reflecting
+  LOG("Enumerate roots.\n");
+
+  memset(root_vectors, 0, 2*positive*rank*sizeof(int));
+  roots_known = 0;
+
+  queue_init(&queue);
+  for(int i = 0; i < rank; i++) {
+    root_vectors[rank*i + i] = 1; // (r_i)_j = delta_ij
+    queue_put(&queue, i);
+    roots_known++;
+  }
+
+  while((current = queue_get(&queue)) != -1) {
+    for(int i = 0; i < rank; i++) {
+      reflect_root_vector(cartan_matrix, rank, i, &root_vectors[rank*current], vector);
+      int j;
+      for(j = 0; j < roots_known; j++)
+	if(compare_root_vectors(rank, &root_vectors[rank*j], vector) == 0)
+	  break;
+      if(j == roots_known) {
+	memcpy(&root_vectors[rank*roots_known], vector, rank*sizeof(int));
+	queue_put(&queue, roots_known);
+	roots_known++;
+      }
+    }
+  }
+
+  ERROR(roots_known != 2*positive, "Number of roots does not match!\n");
+
+  // sort roots and restrict to positives
+  LOG("Sort roots.\n");
+
+  qsort_r(root_vectors, 2*positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
+  memcpy(root_vectors, &root_vectors[positive*rank], positive*rank*sizeof(int)); // this just copies the second part of the list onto the first; source and destination are disjoint!
+
+  // generate root_mapping, which gives the action of root reflections on positive roots (-1 if result is not a positive root)
+  LOG("Compute root reflections.\n");
+
+  for(int i = 0; i < positive; i++) {
+    for(int j = 0; j < rank; j++) {
+      reflect_root_vector(cartan_matrix, rank, j, &root_vectors[rank*i], vector);
+      root_mapping[i*rank+j] =
+	search(vector, root_vectors, positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
+    }
+  }
+
+  // find simple roots in the list
+  LOG("Find simple roots.\n");
+
+  for(int i = 0; i < rank; i++) {
+    memset(vector, 0, rank*sizeof(int));
+    vector[i] = 1;
+    simple_roots[i] = search(vector, root_vectors, positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
+  }
+
+  // enumerate weyl group elements using difference sets
+  LOG("Enumerate Weyl group elements.\n");
+
+  nextids[0] = 0;
+  nextids_count = 1;
+  elements = 0;
+  for(int len = 0; len <= positive; len++) {
+    length_elements = 0;
+
+    // find unique ids in edges added in the last iteration
+    qsort(nextids, nextids_count, sizeof(weylid_t), compare_weylid);
+    for(int i = 0; i < nextids_count; i++)
+      if(i == 0 || nextids[i] != nextids[i-1])
+	ids[elements + length_elements++] = nextids[i];
+
+    // add new edges
+    nextids_count = 0;
+    for(int i = elements; i < elements + length_elements; i++)
+      for(int j = 0; j < rank; j++) {
+	edges[i*rank+j] = multiply_generator(j, ids[i], simple_roots, root_mapping, rank, positive);
+	if(!(ids[i] & BIT(simple_roots[j]))) // the new element is longer then the old one
+	  nextids[nextids_count++] = edges[i*rank+j];
+      }
+
+    elements += length_elements;
+  }
+
+  // translate the ids to list positions (i.e. local continuous ids)
+  LOG("Reorder Weyl group elements.\n");
+
+  for(int i = 0; i < order; i++) {
+    lookup[i].id = ids[i];
+    lookup[i].position = i;
+  }
+  qsort(lookup, order, sizeof(weylid_lookup_t), compare_weylid_lookup);
+
+  // fill in results
+  LOG("Compute left multiplication.\n");
+
+  for(int i = 0; i < order; i++) {
+    group[i].id = ids[i];
+    for(int j = 0; j < rank; j++)
+      group[i].left[j] = group + lookup_id(edges[i*rank+j], lookup, order);
+  }
+
+  // free temporary stuff
+
+  free(cartan_matrix);
+  free(root_vectors);
+  free(vector);
+  free(root_mapping);
+  free(simple_roots);
+  free(ids);
+  free(edges);
+  free(nextids);
+  free(lookup);
+}
 
 // glibc search function, but with user pointer and returning index (or -1 if not found)
 static int search (const void *key, const void *base, size_t nmemb, size_t size, int (*compar) (const void *, const void *, void *), void *arg)
@@ -176,46 +328,46 @@ int weyl_order(semisimple_type_t type)
   return order;
 }
 
-int weyl_hyperplanes(semisimple_type_t type)
+int weyl_positive(semisimple_type_t type)
 {
-  int hyperplanes = 0;
+  int positive = 0;
 
   for(int i = 0; i < type.n; i++) {
     ERROR(!weyl_exists(type.factors[i]), "A Weyl group of type %c%d does not exist!\n", type.factors[i].series, type.factors[i].rank);
 
     switch(type.factors[i].series) {
     case 'A':
-      hyperplanes += (type.factors[i].rank * (type.factors[i].rank + 1)) / 2;
+      positive += (type.factors[i].rank * (type.factors[i].rank + 1)) / 2;
       break;
 
     case 'B': case 'C':
-      hyperplanes += type.factors[i].rank * type.factors[i].rank;
+      positive += type.factors[i].rank * type.factors[i].rank;
       break;
 
     case 'D':
-      hyperplanes += type.factors[i].rank * (type.factors[i].rank - 1);
+      positive += type.factors[i].rank * (type.factors[i].rank - 1);
       break;
 
     case 'E':
       if(type.factors[i].rank == 6)
-	hyperplanes += 36;
+	positive += 36;
       else if(type.factors[i].rank == 7)
-	hyperplanes += 63;
+	positive += 63;
       else if(type.factors[i].rank == 8)
-	hyperplanes += 120;
+	positive += 120;
       break;
 
     case 'F':
-      hyperplanes += 24;
+      positive += 24;
       break;
 
     case 'G':
-      hyperplanes += 6;
+      positive += 6;
       break;
     }
   }
 
-  return hyperplanes;
+  return positive;
 }
 
 int weyl_opposition(semisimple_type_t type, int simple_root)
@@ -313,157 +465,354 @@ void weyl_cartan_matrix(semisimple_type_t type, int *m)
   free(A);
 }
 
-/************ memory allocation ********************/
+/************ weyl_generate etc. ********************/
 
-weylgroup_element_t *weyl_alloc(semisimple_type_t type)
+static weylgroup_element_t* apply_word(int *word, int len, weylgroup_element_t *current)
+{
+  for(int k = len - 1; k >= 0; k--) // apply group element from right to left
+    current = current->left[word[k]];
+
+  return current;
+}
+
+static weylgroup_element_t* apply_word_reverse(int *word, int len, weylgroup_element_t *current)
+{
+  for(int k = 0; k < len; k++) // apply group element from left to right (i.e. apply inverse)
+    current = current->left[word[k]];
+
+  return current;
+}
+
+weylgroup_t *weyl_generate(semisimple_type_t type)
 {
   int rank = weyl_rank(type);
   int order = weyl_order(type);
-
-  int *left = (int*)malloc(rank*order*sizeof(int));
-  int *right = (int*)malloc(rank*order*sizeof(int));
-  weylgroup_element_t *group = (weylgroup_element_t*)malloc(order*sizeof(weylgroup_element_t));
-
-  for(int i = 0; i < order; i++) {
-    group[i].left = &left[i*rank];
-    group[i].right = &right[i*rank];
-  }
-
-  return group;
-}
-
-void weyl_free(weylgroup_element_t *x)
-{
-  free(x[0].left);
-  free(x[0].right);
-  free(x);
-}
-
-void weyl_generate(semisimple_type_t type, weylgroup_element_t *group)
-{
-  int rank, order, positive;
-  queue_t queue;
-  int current;
-  int roots_known, elements, length_elements, nextids_count;
-  int *cartan_matrix;
-  int *root_vectors;
-  int *vector;
-  int *simple_roots;
-  int *root_mapping;
-  weylid_t *ids, *edges, *nextids;
-  weylid_lookup_t *lookup;
-
-  rank = weyl_rank(type);
-  order = weyl_order(type);
-  positive = weyl_hyperplanes(type);
+  int positive = weyl_positive(type);
 
   ERROR(positive > 64, "We can't handle root systems with more than 64 positive roots!\n");
 
-  cartan_matrix =      (int*)malloc(rank*rank      *sizeof(int));
-  root_vectors =       (int*)malloc(2*positive*rank*sizeof(int));
-  vector =             (int*)malloc(rank           *sizeof(int));
-  root_mapping =       (int*)malloc(positive*rank  *sizeof(int));
-  simple_roots =       (int*)malloc(rank           *sizeof(int));
-  ids =           (weylid_t*)malloc(order          *sizeof(weylid_t));
-  edges =         (weylid_t*)malloc(rank*order     *sizeof(weylid_t));
-  nextids =       (weylid_t*)malloc(rank*order     *sizeof(weylid_t));
-  lookup = (weylid_lookup_t*)malloc(order          *sizeof(weylid_lookup_t));
+  // allocate result
 
-  weyl_cartan_matrix(type, cartan_matrix);
+  weylgroup_element_t *group = (weylgroup_element_t*)malloc(order*sizeof(weylgroup_element_t));
+  weylgroup_t *result = malloc(sizeof(weylgroup_t));
+  result->type = type;
+  result->elements = group;
+  result->lists = (weylgroup_element_t**)malloc(2*order*rank*sizeof(weylgroup_element_t*));
 
-  // enumerate roots
-  memset(root_vectors, 0, 2*positive*rank*sizeof(int));
-
-  // first the simple roots
-  queue_init(&queue);
-  for(int i = 0; i < rank; i++) {
-    root_vectors[rank*i + i] = 1;
-    queue_put(&queue, i);
+  for(int i = 0; i < order; i++) {
+    group[i].left = result->lists + 2*i*rank;
+    group[i].right = result->lists + (2*i+1)*rank;
+    group[i].coset = (doublecoset_t*)0;
   }
 
-  // and then we get all others by reflecting
-  roots_known = rank;
-  while((current = queue_get(&queue)) != -1) {
-    for(int i = 0; i < rank; i++) {
-      reflect_root_vector(cartan_matrix, rank, i, &root_vectors[rank*current], vector);
-      int j;
-      for(j = 0; j < roots_known; j++)
-	if(compare_root_vectors(rank, &root_vectors[rank*j], vector) == 0)
-	  break;
-      if(j == roots_known) {
-	memcpy(&root_vectors[rank*roots_known], vector, rank*sizeof(int));
-	queue_put(&queue, roots_known);
-	roots_known++;
-      }
-    }
+  // the main part
+  LOG("Start generating Weyl group.\n");
+
+  generate_left_and_ids(type, group);
+
+  // word length is just the number of 1s in the binary id
+  LOG("Find word lengths.\n");
+
+  for(int i = 0; i < order; i++) {
+    group[i].wordlength = 0;
+    for(int j = 0; j < positive; j++)
+      if(group[i].id & BIT(j))
+	group[i].wordlength++;
   }
 
-  ERROR(roots_known != 2*positive, "Number of roots does not match!\n")
+  // allocate letters
 
-  // sort roots and restrict to positives
-  qsort_r(root_vectors, 2*positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
-  memcpy(root_vectors, &root_vectors[positive*rank], positive*rank*sizeof(int));
+  int total_wordlength = 0;
+  for(int i = 0; i < order; i++)
+    total_wordlength += group[i].wordlength;
+  result->letters = (int*)malloc(total_wordlength*sizeof(int));
+  total_wordlength = 0;
+  for(int i = 0; i < order; i++) {
+    group[i].word = result->letters + total_wordlength;
+    total_wordlength += group[i].wordlength;
+  }
 
-  for(int i = 0; i < positive; i++) {
+  // find shortest words (using that the elements are already ordered by word length)
+  LOG("Find shortest words.\n");
+
+  memset(result->letters, -1, total_wordlength*sizeof(int));
+  for(int i = 0; i < order - 1; i++) {
+    weylgroup_element_t *this = &group[i];
     for(int j = 0; j < rank; j++) {
-      reflect_root_vector(cartan_matrix, rank, j, &root_vectors[rank*i], vector);
-      root_mapping[i*rank+j] =
-	search(vector, root_vectors, positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
+      weylgroup_element_t *that = group[i].left[j];
+      if(that->wordlength > this->wordlength && that->word[0] == -1) {
+	memcpy(that->word + 1, this->word, this->wordlength*sizeof(int));
+	that->word[0] = j;
+      }
     }
   }
 
-  // where in the list are the simple roots?
+  // generate right edges
+  LOG("Compute right multiplication.\n");
+
+  for(int i = 0; i < order; i++)
+    for(int j = 0; j < rank; j++)
+      group[i].right[j] = apply_word(group[i].word, group[i].wordlength, group[0].left[j]);
+
+  // find opposites
+  LOG("Find opposites.\n");
+
+  weylgroup_element_t *longest = &group[order-1];
+  for(int i = 0; i < order; i++)
+    group[i].opposite = apply_word(longest->word, longest->wordlength, &group[i]);
+
+  // check for root reflections
+  LOG("Find root reflections.\n");
+
+  for(int i = 0; i < order; i++)
+    group[i].is_root_reflection = 0;
+  for(int i = 0; i < order; i++)
+    for(int j = 0; j < rank; j++) // we want to calculate word^{-1} * j * word; this is a root reflection
+      apply_word_reverse(group[i].word, group[i].wordlength, group[i].left[j]) -> is_root_reflection = 1;
+
+  return result;
+}
+
+void weyl_destroy(weylgroup_t *group)
+{
+  free(group->elements);
+  free(group->lists);
+  free(group->letters);
+  free(group);
+}
+
+doublequotient_t *weyl_generate_bruhat(semisimple_type_t type, int left_invariance, int right_invariance)
+{
+  int rank = weyl_rank(type);
+  int order = weyl_order(type);
+  int positive = weyl_positive(type);
+  int count;
+
+  int is_minimum, is_maximum;
+
+  weylgroup_t *wgroup = weyl_generate(type);
+  weylgroup_element_t *group = wgroup->elements;
+  doublecoset_t *cosets;
+
   for(int i = 0; i < rank; i++) {
-    memset(vector, 0, rank*sizeof(int));
-    vector[i] = 1;
-    simple_roots[i] = search(vector, root_vectors, positive, rank*sizeof(int), compare_root_vectors_qsort, &rank);
+    int oppi = weyl_opposition(type, i);
+    if(left_invariance & BIT(i) && !(left_invariance & BIT(oppi)) ||
+       left_invariance & BIT(oppi) && !(left_invariance & BIT(i)))
+      ERROR(1, "The specified left invariance is not invariant under the opposition involution!\n");
   }
 
-  // enumerate weyl group elements using difference sets
-  nextids[0] = 0;
-  nextids_count = 1;
-  elements = 0;
-  for(int len = 0; len <= positive; len++) {
-    length_elements = 0;
+  doublequotient_t *result = (doublequotient_t*)malloc(sizeof(doublequotient_t));
+  result->type = type;
+  result->left_invariance = left_invariance;
+  result->right_invariance = right_invariance;
+  result->group = wgroup->elements;
+  result->grouplists = wgroup->lists;
+  result->groupletters = wgroup->letters;
 
-    // find unique ids in edges added in the last iteration
-    qsort(nextids, nextids_count, sizeof(weylid_t), compare_weylid);
-    for(int i = 0; i < nextids_count; i++)
-      if(i == 0 || nextids[i] != nextids[i-1])
-	ids[elements + length_elements++] = nextids[i];
+  free(wgroup); // dissolved in result and not needed anymore
 
-    // add new edges
-    nextids_count = 0;
-    for(int i = elements; i < elements + length_elements; i++)
-      for(int j = 0; j < rank; j++) {
-	edges[i*rank+j] = multiply_generator(j, ids[i], simple_roots, root_mapping, rank, positive);
-	if(!(ids[i] & BIT(simple_roots[j]))) // the new element is longer then the old one
-	  nextids[nextids_count++] = edges[i*rank+j];
+  // count cosets by finding the minimum length element in every coset
+  LOG("Count cosets.\n");
+
+  count = 0;
+  for(int i = 0; i < order; i++) {
+    is_minimum = 1;
+    for(int j = 0; j < rank; j++)
+      if(left_invariance  & BIT(j) && group[i].left[j]->wordlength  < group[i].wordlength ||
+	 right_invariance & BIT(j) && group[i].right[j]->wordlength < group[i].wordlength)
+	is_minimum = 0;
+    if(is_minimum)
+      count++;
+  }
+  result->count = count;
+
+  // alloc more stuff
+
+  cosets = result->cosets = (doublecoset_t*)malloc(count*sizeof(doublecoset_t));
+  for(int i = 0; i < count; i++) {
+    cosets[i].bruhat_lower = cosets[i].bruhat_higher = (doublecoset_list_t*)0;
+  }
+  result->lists = (doublecoset_list_t*)malloc(2*count*positive*sizeof(doublecoset_list_t)); // 2 times, for bruhat lower and higher
+
+  // find minima (basically same code as above)
+  LOG("Find minimal length elements in cosets.\n");
+
+  count = 0;
+  for(int i = 0; i < order; i++) {
+    is_minimum = 1;
+    for(int j = 0; j < rank; j++)
+      if(left_invariance  & BIT(j) && group[i].left[j]->wordlength  < group[i].wordlength ||
+	 right_invariance & BIT(j) && group[i].right[j]->wordlength < group[i].wordlength)
+	is_minimum = 0;
+    if(is_minimum) {
+      cosets[count].min = &group[i];
+      group[i].coset = &cosets[count];
+      count++;
+    }
+  }
+
+  // generate quotient map
+  LOG("Generate quotient map.\n");
+
+  for(int i = 0; i < order; i++) {
+    for(int j = 0; j < rank; j++) {
+      if(left_invariance & BIT(j) && group[i].left[j]->wordlength > group[i].wordlength)
+	group[i].left[j]->coset = group[i].coset;
+      if(right_invariance & BIT(j) && group[i].right[j]->wordlength > group[i].wordlength)
+	group[i].right[j]->coset = group[i].coset;
+    }
+  }
+
+  // find maxima
+  LOG("Find maximal length elements.\n");
+
+  for(int i = 0; i < order; i++) {
+    is_maximum = 1;
+    for(int j = 0; j < rank; j++)
+      if(left_invariance  & BIT(j) && group[i].left[j]->wordlength  > group[i].wordlength ||
+	 right_invariance & BIT(j) && group[i].right[j]->wordlength > group[i].wordlength)
+	is_maximum = 0;
+    if(is_maximum) {
+      group[i].coset->max = &group[i];
+    }
+  }
+
+  // opposites
+  LOG("Find opposites.\n");
+
+  for(int i = 0; i < count; i++)
+    cosets[i].opposite = cosets[i].min->opposite->coset;
+
+  // bruhat order
+  LOG("Find bruhat order.\n");
+
+  int edgecount = 0;
+  for(int i = 0; i < order; i++) {
+    if(group[i].is_root_reflection) {
+      for(int j = 0; j < count; j++) {
+	weylgroup_element_t *this = cosets[j].min;
+	weylgroup_element_t *that = apply_word(group[i].word, group[i].wordlength, cosets[j].min);
+	if(this->wordlength > that->wordlength) { // this is higher in bruhat order than that
+	  doublecoset_list_t *new = &result->lists[edgecount++];
+	  new->next = this->coset->bruhat_lower;
+	  this->coset->bruhat_lower = new;
+	  new->to = that->coset;
+	}
+      }
+    }
+  }
+
+  // transitive reduction
+  LOG("Perform transitive reduction.\n");
+
+  doublecoset_t *offset = &cosets[0];
+  doublecoset_t *origin;
+  doublecoset_list_t *current;
+  doublecoset_list_t *prev;
+  queue_t queue;
+  int cur;
+  int *seen = malloc(count*sizeof(int));
+  for(int i = 0; i < count; i++) {
+    memset(seen, 0, count*sizeof(int));
+    queue_init(&queue);
+
+    for(int len = 1; len <= cosets[i].min->wordlength; len++) {
+
+      // remove all edges originating from i of length len which connect to something already seen using shorter edges
+      origin = &cosets[i];
+      prev = (doublecoset_list_t*)0;
+
+      for(current = origin->bruhat_lower; current; current = current->next) {
+	if(origin->min->wordlength - current->to->min->wordlength != len) {
+	  prev = current;
+	} else if(seen[current->to - offset]) {
+	  if(prev)
+	    prev->next = current->next;
+	  else
+	    origin->bruhat_lower = current->next;
+	} else {
+	  prev = current;
+	  seen[current->to - offset] = 1;
+	  queue_put(&queue, current->to - offset);
+	}
       }
 
-    elements += length_elements;
+      // see which nodes we can reach using only edges up to length len, mark them as seen
+      while((cur = queue_get(&queue)) != -1) {
+	current = (cur + offset)->bruhat_lower;
+	for(current = (cur+offset)->bruhat_lower; current; current = current->next) {
+	  if(!seen[current->to - offset]) {
+	    seen[current->to - offset] = 1;
+	    queue_put(&queue, current->to - offset);
+	  }
+	}
+      }
+    }
   }
 
-  // translate the ids to list positions (i.e. local continuous ids)
-  for(int i = 0; i < order; i++) {
-    lookup[i].id = ids[i];
-    lookup[i].position = i;
-  }
-  qsort(lookup, order, sizeof(weylid_lookup_t), compare_weylid_lookup);
+  free(seen);
 
-  for(int i = 0; i < order; i++) {
-    group[i].id = ids[i];
-    for(int j = 0; j < rank; j++)
-      group[i].left[j] = lookup_id(edges[i*rank+j], lookup, order);
+  // reverse bruhat order
+  LOG("Revert bruhat order.\n");
+
+  for(int i = 0; i < count; i++) {
+    for(current = cosets[i].bruhat_lower; current; current = current->next) {
+      doublecoset_list_t *new = &result->lists[edgecount++];
+      new->to = &cosets[i];
+      new->next = current->to->bruhat_higher;
+      current->to->bruhat_higher = new;
+    }
   }
 
-  free(cartan_matrix);
-  free(root_vectors);
-  free(vector);
-  free(root_mapping);
-  free(simple_roots);
-  free(ids);
-  free(edges);
-  free(nextids);
-  free(lookup);
+  // sort opposites and rewrite everything
+  LOG("Sort opposites.\n");
+
+  int *old2newindices = (int*)malloc(count*sizeof(int));
+  int *new2oldindices = (int*)malloc(count*sizeof(int));
+  doublecoset_t *oldcosets = (doublecoset_t*)malloc(count*sizeof(doublecoset_t));
+  memcpy(oldcosets, cosets, count*sizeof(doublecoset_t));
+
+  int j = 0;
+  for(int i = 0; i < count; i++)
+    if(&cosets[i] < cosets[i].opposite) {
+      old2newindices[i] = j;
+      old2newindices[cosets[i].opposite - cosets] = count-1-j;
+      j++;
+    }
+  for(int i = 0; i < count; i++)
+    if(i == cosets[i].opposite - cosets)
+      old2newindices[i] = j++;
+
+  for(int i = 0; i < count; i++)
+    new2oldindices[old2newindices[i]] = i;
+
+  for(int i = 0; i < count; i++) {
+    cosets[i].min = oldcosets[new2oldindices[i]].min;
+    cosets[i].max = oldcosets[new2oldindices[i]].max;
+    cosets[i].opposite = old2newindices[oldcosets[new2oldindices[i]].opposite - cosets] + cosets;
+    cosets[i].bruhat_lower = oldcosets[new2oldindices[i]].bruhat_lower;
+    cosets[i].bruhat_higher = oldcosets[new2oldindices[i]].bruhat_higher;
+    for(current = cosets[i].bruhat_lower; current; current = current -> next)
+      current->to = old2newindices[current->to - cosets] + cosets;
+    for(current = cosets[i].bruhat_higher; current; current = current -> next)
+      current->to = old2newindices[current->to - cosets] + cosets;
+  }
+  for(int i = 0; i < order; i++)
+    group[i].coset = old2newindices[group[i].coset - cosets] + cosets;
+
+  free(old2newindices);
+  free(new2oldindices);
+  free(oldcosets);
+
+  return result;
+}
+
+void weyl_destroy_bruhat(doublequotient_t *dq)
+{
+  free(dq->group);
+  free(dq->grouplists);
+  free(dq->groupletters);
+  free(dq->cosets);
+  free(dq->lists);
+  free(dq);
 }
diff --git a/weyl.h b/weyl.h
index fba78cf..55b8442 100644
--- a/weyl.h
+++ b/weyl.h
@@ -3,34 +3,92 @@
 
 #include <inttypes.h>
 
-typedef struct {
-  char series;
-  int rank;
-} simple_type_t;
-
-typedef struct {
-  int n;
-  simple_type_t *factors;
-} semisimple_type_t;
+struct _simple_type;
+struct _semisimple_type;
+struct _weylgroup_element;
+struct _weylgroup;
+struct _doublecoset;
+struct _doublecoset_list;
+struct _doublequotient;
 
 typedef uint64_t weylid_t;
+typedef struct _simple_type simple_type_t;
+typedef struct _semisimple_type semisimple_type_t;
+typedef struct _weylgroup_element weylgroup_element_t;
+typedef struct _weylgroup weylgroup_t;
+typedef struct _doublecoset doublecoset_t;
+typedef struct _doublecoset_list doublecoset_list_t;
+typedef struct _doublequotient doublequotient_t;
 
-typedef struct {
-  int *left;
-  int *right;
-  int opposite;
+/***************************** structures *******************************/
+
+struct _simple_type {
+  char series;
+  int rank;
+};
+
+struct _semisimple_type {
+  int n;
+  simple_type_t *factors;
+};
+
+struct _weylgroup_element {
+  int *word;
+  int wordlength;
+  weylgroup_element_t **left;
+  weylgroup_element_t **right;
+  weylgroup_element_t *opposite;
+  int is_root_reflection; // boolean value
   weylid_t id;
-} weylgroup_element_t;
+
+  // only set if quotient is generated
+  doublecoset_t *coset;
+};
+
+struct _weylgroup {
+  semisimple_type_t type;
+  weylgroup_element_t *elements;
+  weylgroup_element_t **lists;
+  int *letters;
+};
+
+struct _doublecoset {
+  doublecoset_list_t *bruhat_lower;
+  doublecoset_list_t *bruhat_higher;
+  doublecoset_t *opposite;
+  weylgroup_element_t *max;
+  weylgroup_element_t *min;
+};
+
+struct _doublecoset_list {
+  doublecoset_t *to;
+  doublecoset_list_t *next;
+};
+
+struct _doublequotient {
+  semisimple_type_t type;
+  int left_invariance;                // bitmask with rank bits
+  int right_invariance;
+  int count;                          // number of cosets
+  doublecoset_t *cosets;
+  weylgroup_element_t *group;
+  doublecoset_list_t *lists;          // only for memory allocation / freeing
+  weylgroup_element_t **grouplists;   // only for memory allocation / freeing
+  int *groupletters;                  // only for memory allocation / freeing
+};
+
+/***************************** functions **************************************/
 
 int weyl_rank(semisimple_type_t type);
 int weyl_order(semisimple_type_t type);
-int weyl_hyperplanes(semisimple_type_t type);
+int weyl_positive(semisimple_type_t type);
 void weyl_cartan_matrix(semisimple_type_t type, int *m);
 int weyl_opposition(semisimple_type_t type, int simple_root);
 
-weylgroup_element_t *weyl_alloc(semisimple_type_t type);
-void weyl_free(weylgroup_element_t *x);
+weylgroup_t *weyl_generate(semisimple_type_t type);
+void weyl_destroy(weylgroup_t *group);
 
-void weyl_generate(semisimple_type_t type, weylgroup_element_t *group);
+doublequotient_t *weyl_generate_bruhat(semisimple_type_t type, int left_invariance, int right_invariance);
+void weyl_destroy_bruhat(doublequotient_t *dq);
 
 #endif