New Weyl group algorithm

Makefile

@@ -1,15 +1,18 @@
-HEADERS=coxeter.h thickenings.h queue.h bitvec.h
-OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99 -D_GNU_SOURCE -Winline
-#OPTIONS=-m64 -march=native -O0 -g -std=gnu99
-#OPTIONS=-O3 -m64 -march=native -funroll-loops -fno-inline -std=gnu99 -pg
+HEADERS=weyl.h thickenings.h queue.h bitvec.h
+
+SPECIAL_OPTIONS=-O0 -g
+#SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
+#SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
+
+OPTIONS=-m64 -march=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
 
 all: generate process
 
-generate: generate.o coxeter.o thickenings.o
-	gcc $(OPTIONS) -o generate generate.o thickenings.o coxeter.o -lgsl -lcblas
+generate: generate.o weyl.o thickenings.o
+	gcc $(OPTIONS) -o generate generate.o thickenings.o weyl.o
 
-process: process.o coxeter.o thickenings.o
-	gcc $(OPTIONS) -o process process.o thickenings.o coxeter.o -lgsl -lcblas
+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
@@ -20,8 +23,8 @@ process.o: process.c $(HEADERS)
 thickenings.o: thickenings.c $(HEADERS)
 	gcc $(OPTIONS) -c thickenings.c
 
-coxeter.o: coxeter.c $(HEADERS)
-	gcc $(OPTIONS) -c coxeter.c
+weyl.o: weyl.c $(HEADERS)
+	gcc $(OPTIONS) -c weyl.c
 
 clean:
-	rm -f generate process thickenings.o coxeter.o generate.o process.o
+	rm -f generate process thickenings.o weyl.o generate.o process.o

generate.c

@@ -1,5 +1,5 @@
 #include "thickenings.h"
-#include "coxeter.h"
+#include "weyl.h"
 #include "queue.h"
 
 #include <strings.h>
@@ -64,13 +64,13 @@ int main(int argc, const char *argv[])
   graph = graph_alloc(type);
   cosets = prepare_simplified_graph(type, left_invariance, right_invariance, graph);
 
-  ERROR(cosets < 0, "The left invariance is not preserved by the opposition involution: %d %d!\n", left_invariance, opposition_involution(type, left_invariance));
+  ERROR(cosets < 0, "The left invariance is not preserved by the opposition involution!\n");
 
   // print stuff
 
-  rank = coxeter_rank(type);                // number of simple roots
-  order = coxeter_order(type);              // number of Weyl group elements
-  hyperplanes = coxeter_hyperplanes(type);  // number of positive roots
+  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");

process.c

@@ -3,7 +3,7 @@
 #include <sys/stat.h>
 
 #include "thickenings.h"
-#include "coxeter.h"
+#include "weyl.h"
 #include "queue.h"
 
 int main(int argc, const char *argv[])
@@ -35,7 +35,7 @@ int main(int argc, const char *argv[])
   fread(&left_invariance, sizeof(simple_type_t), type.n, infile);
   fread(&right_invariance, sizeof(simple_type_t), type.n, infile);
 
-  rank = coxeter_rank(type);
+  rank = weyl_rank(type);
   graph = graph_alloc(type);
   cosets = prepare_simplified_graph(type, left_invariance, right_invariance, graph);
 

thickenings.c

@@ -5,10 +5,9 @@
 #include <memory.h>
 
 #include "thickenings.h"
-#include "coxeter.h"
+#include "weyl.h"
 #include "queue.h"
 
-
 char *alphabetize(int *word, int len, const char *alphabet, char *buffer)
 {
   if(len == 0) {
@@ -66,20 +65,20 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
   int edgelist_count, hyperplane_count;
   int current;
 
-  int *graph_data;
+  weylgroup_element_t *graph_data;
   node_t *graph_unsorted;
   int *ordering, *reverse_ordering, *seen;
 
   // initialize
 
-  rank = coxeter_rank(type);
-  order = coxeter_order(type);
-  hyperplanes = coxeter_hyperplanes(type);
+  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 = (int*)malloc(order*rank*sizeof(int));
+  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));
@@ -94,11 +93,15 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
 
   // get coxeter graph
 
-  generate_coxeter_graph(type, graph_data);
+  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*rank];
+    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
 
@@ -115,6 +118,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Wordlengths calculated.\n");
+
   // sort by wordlength
 
   for(int i = 0; i < order; i++)
@@ -123,12 +128,15 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
   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 and wordlength so far, so just copy these
+    // 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
   }
 
+  fprintf(stderr, "Sorted by wordlength.\n");
+
   // find words
 
   for(int i = 0; i < order; i++)
@@ -146,6 +154,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Shortest words found.\n");
+
   // generate right edges
 
   for(int i = 0; i < order; i++) {
@@ -158,6 +168,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Right edges generated.\n");
+
   // find opposites
 
   node_t *longest = &graph[order-1];
@@ -168,6 +180,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     graph[i].opposite = current;
   }
 
+  fprintf(stderr, "Opposites found.\n");
+
   // enumerate hyperplanes
 
   hyperplane_count = 0;
@@ -189,6 +203,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Hyperplanes enumerated.\n");
+
   // generate folding order
 
   edgelist_count = 0;
@@ -213,6 +229,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Bruhat order generated.\n");
+
   // remove redundant edges
 
   for(int i = 0; i < order; i++) {
@@ -254,6 +272,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Redundant edges removed.\n");
+
   // reverse folding order
 
   edgelist_count = 0;
@@ -268,6 +288,8 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
     }
   }
 
+  fprintf(stderr, "Bruhat order reversed.\n");
+
   // additional sorting step to force opposite property (opposite of j is at n - j - 1)
 
   for(int i = 0; i < order; i++)
@@ -295,7 +317,7 @@ void prepare_graph(semisimple_type_t type, node_t *graph)
       edge->to = reverse_ordering[edge->to];
   }
 
-  free(graph_data);
+  weyl_free(graph_data);
   free(graph_unsorted);
   free(ordering);
   free(reverse_ordering);
@@ -332,8 +354,16 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne
   queue_t queue;
   int ncosets;
 
-  if(opposition_involution(type, left) != left)
-    return -1;
+  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];
@@ -343,11 +373,9 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne
   full_graph = graph_alloc(type);
   prepare_graph(type, full_graph);
 
-  // initialize stuff
+  fprintf(stderr, "Full graph generated.\n");
 
-  rank = coxeter_rank(type);
-  order = coxeter_order(type);
-  hyperplanes = coxeter_hyperplanes(type);
+  // initialize stuff
 
   reduced = (int*)malloc(order*sizeof(int));
   group = (int*)malloc(order*sizeof(int));
@@ -478,6 +506,7 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne
   }
 
   // step 8: revert order
+  edgelists_used = 0;
   for(int i = 0; i < ncosets; i++) {
     edge = simplified_graph[i].bruhat_lower;
     while(edge) {
@@ -520,9 +549,9 @@ int prepare_simplified_graph(semisimple_type_t type, unsigned long left, unsigne
 
 node_t *graph_alloc(semisimple_type_t type)
 {
-  int rank = coxeter_rank(type);
-  int order = coxeter_order(type);
-  int hyperplanes = coxeter_hyperplanes(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));
@@ -547,7 +576,7 @@ void graph_free(semisimple_type_t type, node_t *graph)
   free(graph[0].right);
   free(graph[0].word);
 
-  int order = coxeter_order(type);
+  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;
@@ -582,6 +611,7 @@ typedef struct {
   - returns number of balanced ideals found
 
   uses the bitvector functions bv_union, bv_copy, bv_set_range_except, bv_disjoint, bv_next_zero
+
  */
 
 static long enumerate_tree(const enumeration_info_t *info, const bitvec_t *pos, const bitvec_t *neg, int next_neg, int already_known)
@@ -698,6 +728,13 @@ long enumerate_balanced_thickenings(node_t *graph, int size, void (*callback) (c
   }
   free(principal);
 
+  /*
+  for(int i = 0; i < info.size; i++) {
+    fprintf(stderr, "%d: ", i);
+    bv_print_nice(stderr, &info.principal_pos[i], &info.principal_neg[i], -1, info.size/2);
+    fprintf(stderr, "\n");
+    } */
+
   // enumerate balanced ideals
   bitvec_t pos, neg;
   bv_clear(&pos);

thickenings.h

@@ -4,7 +4,7 @@
 #define BV_QWORD_RANK 10
 #include "bitvec.h"
 
-#include "coxeter.h"
+#include "weyl.h"
 
 #define DEBUG(msg, ...) do{fprintf(stderr, msg, ##__VA_ARGS__); }while(0)
 
@@ -16,7 +16,7 @@ typedef struct _edgelist {
   struct _edgelist *next;
 } edgelist_t;
 
-// describes an element of the Coxeter group; only "opposite" and "bruhat_lower" are being used for enumerating thickenings; everything else is just needed for initialization or output
+// 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;
@@ -26,6 +26,7 @@ typedef struct {
   edgelist_t *bruhat_lower;
   edgelist_t *bruhat_higher;
   int is_hyperplane_reflection; // boolean value
+  weylid_t id;
 } node_t;
 
 // printing functions

weyl.c

@@ -0,0 +1,487 @@
+#include "weyl.h"
+#include "queue.h"
+
+#include <stdio.h>
+#include <memory.h>
+#include <stdlib.h>
+
+#define BIT(n) ((uint64_t)1 << (n))
+
+typedef struct {
+  weylid_t id;
+  int position;
+} weylid_lookup_t;
+
+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_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);
+
+/***************** simple helper functions **********************************/
+
+// 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)
+{
+  size_t l, u, idx;
+  const void *p;
+  int comparison;
+
+  l = 0;
+  u = nmemb;
+  while (l < u) {
+    idx = (l + u) / 2;
+    p = (void *) (((const char *) base) + (idx * size));
+    comparison = (*compar) (key, p, arg);
+    if (comparison < 0)
+      u = idx;
+    else if (comparison > 0)
+      l = idx + 1;
+    else
+      return idx;
+  }
+
+  return -1;
+}
+
+// maybe we want a different ordering here?
+static int compare_root_vectors(int rank, const int *x, const int *y)
+{
+  for(int i = 0; i < rank; i++)
+    if(x[i] != y[i])
+      return x[i] - y[i];
+
+  return 0;
+}
+
+static int compare_root_vectors_qsort(const void *x, const void *y, void *arg)
+{
+  return compare_root_vectors(*((int*)arg), x, y);
+}
+
+static int compare_weylid(const void *x, const void *y)
+{
+  weylid_t u = *((weylid_t*)x);
+  weylid_t v = *((weylid_t*)y);
+
+  return u > v ? 1 : u < v ? -1 : 0;
+}
+
+static int compare_weylid_lookup(const void *x, const void *y)
+{
+  weylid_t u = ((weylid_lookup_t*)x)->id;
+  weylid_t v = ((weylid_lookup_t*)y)->id;
+
+  return u > v ? 1 : u < v ? -1 : 0;
+}
+
+static int lookup_id(weylid_t id, weylid_lookup_t *list, int len)
+{
+  weylid_lookup_t key;
+  key.id = id;
+  weylid_lookup_t *p = (weylid_lookup_t*)bsearch(&key, list, len, sizeof(weylid_lookup_t), compare_weylid_lookup);
+  return p->position;
+}
+
+static weylid_t multiply_generator(int s, weylid_t w, const int* simple, const int* mapping, int rank, int positive)
+{
+  weylid_t sw = 0;
+
+  for(int i = 0; i < positive; i++) {
+    if(w & BIT(i))
+      if(mapping[i*rank+s] != -1)
+	sw |= BIT(mapping[i*rank+s]);
+  }
+
+  if(w & BIT(simple[s]))
+    return sw;
+  else
+    return sw | BIT(simple[s]);
+}
+
+static void reflect_root_vector(const int *cartan, int rank, int i, int *old, int *new)
+{
+  memcpy(new, old, rank*sizeof(int));
+  for(int j = 0; j < rank; j++)
+    new[i] -= cartan[i*rank + j]*old[j];
+}
+
+/************* Weyl group infos ************************/
+
+int weyl_rank(semisimple_type_t type)
+{
+  int rank = 0;
+  for(int i = 0; i < type.n; i++)
+    rank += type.factors[i].rank;
+  return rank;
+}
+
+int weyl_order(semisimple_type_t type)
+{
+  int order = 1;
+  for(int i = 0; i < type.n; i++) {
+    switch(type.factors[i].series) {
+
+    case 'A':
+      for(int j = 1; j <= type.factors[i].rank + 1; j++)
+	order *= j;
+      break;
+
+    case 'B': case 'C':
+      for(int j = 1; j <= type.factors[i].rank; j++)
+	order *= 2*j;
+      break;
+
+    case 'D':
+      for(int j = 2; j <= type.factors[i].rank; j++)
+	order *= 2*j;
+      break;
+
+    case 'E':
+      if(type.factors[i].rank == 6)
+	order *= 51840;
+      else if(type.factors[i].rank == 7)
+	order *= 2903040;
+      else if(type.factors[i].rank == 8)
+	order *= 696729600;
+      else
+	ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+      break;
+
+    case 'F':
+      ERROR(type.factors[i].rank != 4, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+      order *= 1152;
+      break;
+
+    case 'G':
+      ERROR(type.factors[i].rank != 2, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+      order *= 12;
+      break;
+
+    default:
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+    }
+  }
+
+  return order;
+}
+
+int weyl_hyperplanes(semisimple_type_t type)
+{
+  int hyperplanes = 0;
+
+  for(int i = 0; i < type.n; i++) {
+    switch(type.factors[i].series) {
+    case 'A':
+      hyperplanes += (type.factors[i].rank * (type.factors[i].rank + 1)) / 2;
+      break;
+
+    case 'B': case 'C':
+      hyperplanes += type.factors[i].rank * type.factors[i].rank;
+      break;
+
+    case 'D':
+      hyperplanes += type.factors[i].rank * (type.factors[i].rank - 1);
+      break;
+
+    case 'E':
+      if(type.factors[i].rank == 6)
+	hyperplanes += 36;
+      else if(type.factors[i].rank == 7)
+	hyperplanes += 63;
+      else if(type.factors[i].rank == 8)
+	hyperplanes += 120;
+      else
+	ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+      break;
+
+    case 'F':
+      ERROR(type.factors[i].rank != 4, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+      hyperplanes += 24;
+      break;
+
+    case 'G':
+      ERROR(type.factors[i].rank != 2, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+      hyperplanes += 6;
+      break;
+
+    default:
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[i].series, type.factors[i].rank);
+    }
+  }
+
+  return hyperplanes;
+}
+
+int weyl_opposition(semisimple_type_t type, int simple_root)
+{
+  int offset = 0;
+  int factor = 0;
+  int r, iota_r;
+
+  for(factor = 0; factor < type.n; factor++)
+    if(simple_root < offset + type.factors[factor].rank)
+      break;
+    else
+      offset += type.factors[factor].rank;
+  r = simple_root - offset;
+
+  switch(type.factors[factor].series) {
+  case 'A':
+    iota_r = type.factors[factor].rank - 1 - r;
+    break;
+
+  case 'B': case 'C':
+    iota_r = r;
+    break;
+
+  case 'D':
+    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
+    break;
+
+  case 'E':
+    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
+    break;
+
+  case 'F':
+    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
+    break;
+
+  case 'G':
+    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
+    break;
+
+  default:
+    ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[factor].series, type.factors[factor].rank);
+  }
+
+  return iota_r + offset;
+}
+
+void weyl_cartan_matrix(semisimple_type_t type, int *m)
+{
+  int offset = 0;
+  int rank = weyl_rank(type);
+
+  int **A = (int**)malloc(rank*sizeof(int*));
+
+  memset(m, 0, rank*rank*sizeof(int));
+  for(int i = 0; i < rank; i++)
+    m[i*rank+i] = 2;
+
+  for(int k = 0; k < type.n; k++) {
+    for(int i = 0; i < type.factors[k].rank; i++)  // A is the submatrix corresponding to the current simple factor
+      A[i] = &m[(i+offset)*rank + offset];
+
+    switch(type.factors[k].series) {
+    case 'A':
+      for(int i = 1; i < type.factors[k].rank; i++) {
+	A[i][i-1] = -1;
+	A[i-1][i] = -1;
+      }
+      break;
+
+    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;
+      }
+      for(int i = 2; i < type.factors[k].rank; i++) {
+	A[i][i-1] = -1;
+	A[i-1][i] = -1;
+      }
+      break;
+
+    case 'C':
+      if(type.factors[k].rank >= 2) {
+	A[0][1] = -2;
+	A[1][0] = -1;
+      }
+      for(int i = 2; i < type.factors[k].rank; i++) {
+	A[i][i-1] = -1;
+	A[i-1][i] = -1;
+      }
+      break;
+
+    case 'D':
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
+      break;
+
+    case 'E':
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
+      break;
+
+    case 'F':
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
+      break;
+
+    case 'G':
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
+      break;
+
+    default:
+      ERROR(1, "A Weyl group of type %c%d does not exist or is not implemented!\n", type.factors[k].series, type.factors[k].rank);
+    }
+
+    offset += type.factors[k].rank;
+  }
+
+  free(A);
+}
+
+/************ memory allocation ********************/
+
+weylgroup_element_t *weyl_alloc(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);
+
+  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));
+
+  weyl_cartan_matrix(type, cartan_matrix);
+
+  // 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);
+  }
+
+  // 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++;
+      }
+    }
+  }
+
+  ERROR(roots_known != 2*positive, "Number of roots does not match!\n")
+
+  // 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));
+
+  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);
+    }
+  }
+
+  // where in the list are the simple roots?
+  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
+  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)
+  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);
+
+  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);
+  }
+
+  free(cartan_matrix);
+  free(root_vectors);
+  free(vector);
+  free(root_mapping);
+  free(simple_roots);
+  free(ids);
+  free(edges);
+  free(nextids);
+  free(lookup);
+}

weyl.h

@@ -0,0 +1,36 @@
+#ifndef WEYL_H
+#define WEYL_H
+
+#include <inttypes.h>
+
+typedef struct {
+  char series;
+  int rank;
+} simple_type_t;
+
+typedef struct {
+  int n;
+  simple_type_t *factors;
+} semisimple_type_t;
+
+typedef uint64_t weylid_t;
+
+typedef struct {
+  int *left;
+  int *right;
+  int opposite;
+  weylid_t id;
+} weylgroup_element_t;
+
+int weyl_rank(semisimple_type_t type);
+int weyl_order(semisimple_type_t type);
+int weyl_hyperplanes(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);
+
+void weyl_generate(semisimple_type_t type, weylgroup_element_t *group);
+
+#endif