copy .gitignore and remove unnecessary mtune option

.gitignore

@@ -0,0 +1,6 @@
+*.o
+enumerate
+idealbounds
+graph
+D2n
+dominant_weights

D2n.c

@@ -0,0 +1,133 @@
+#include "thickenings.h"
+#include "weyl.h"
+#include "queue.h"
+
+#include <strings.h>
+#include <stdio.h>
+#include <inttypes.h>
+
+int main(int argc, const char *argv[])
+{
+	if(argc < 2) {
+		fprintf(stderr, "Rank argument required.\n");
+		return 1;
+	}
+	int rank = atoi(argv[1]);
+	if(rank <= 0 || rank > 1000) {
+		fprintf(stderr, "Rank must be a small positive integer.\n");
+		return 1;
+	}
+	if(rank % 2) {
+		fprintf(stderr, "Rank must be even.\n"); // opposition involution not trivial otherwise
+		return 1;
+	}
+
+	doublequotient_t *dq = (doublequotient_t*)malloc(sizeof(doublequotient_t));
+	simple_type_t type;
+	type.series = 'D';
+	type.rank = rank;
+	dq->type.n = 1;
+	dq->type.factors = &type;
+	dq->left_invariance = (1<<rank) - 2;
+	dq->right_invariance = 0;
+	dq->group = 0;
+	dq->grouplists = 0;
+	dq->groupletters = 0;
+	dq->count = 1 << (rank - 1);
+	dq->cosets = (doublecoset_t*)malloc(dq->count*sizeof(doublecoset_t));
+	dq->lists = (doublecoset_list_t*)malloc(2*dq->count*rank*sizeof(doublecoset_list_t)); // conservative estimate
+	int nlists = 0;
+
+	int *bitmask = malloc((1<<(rank-1)) * sizeof(int));
+	int *bitmask_reverse = malloc((1<<rank) * sizeof(int));
+	int index = 0;
+
+	LOG("Prepare.\n");
+
+	for(int i = 0; i < dq->count; i++) {
+		dq->cosets[i].index = i;
+		dq->cosets[i].bruhat_lower = (doublecoset_list_t*)0;
+		dq->cosets[i].bruhat_higher = (doublecoset_list_t*)0;
+		dq->cosets[i].min = (weylgroup_element_t*)0;
+		dq->cosets[i].max = (weylgroup_element_t*)0;
+	}
+
+	LOG("Create bitmasks.\n");
+
+	for(uint64_t i = 0; i < (1 << rank); i++) {
+		if(__builtin_popcountll(i) % 2 == 1)
+			bitmask_reverse[i] = -1;
+		else {
+			bitmask[index] = i;
+			bitmask_reverse[i] = index++;
+		}
+	}
+
+	LOG("Generate bruhat order.\n");
+
+	for(int i = 0; i < dq->count; i++) {
+		for(int j = 0; j < rank - 1; j++) {
+			if(!(bitmask[i] & BIT(j)) && (bitmask[i] & BIT(j+1))) {
+				int lowerind = bitmask_reverse[bitmask[i] ^ (BIT(j) | BIT(j+1))];
+				dq->lists[nlists].next = dq->cosets[i].bruhat_lower;
+				dq->cosets[i].bruhat_lower = &dq->lists[nlists];
+				dq->cosets[i].bruhat_lower->to = &dq->cosets[lowerind];
+				nlists++;
+			}
+		}
+		if((bitmask[i] & BIT(0)) && (bitmask[i] & BIT(1))) {
+			int lowerind = bitmask_reverse[bitmask[i] & ~(BIT(0) | BIT(1))];
+			dq->lists[nlists].next = dq->cosets[i].bruhat_lower;
+			dq->cosets[i].bruhat_lower = &dq->lists[nlists];
+			dq->cosets[i].bruhat_lower->to = &dq->cosets[lowerind];
+			nlists++;
+		}
+	}
+
+	LOG("Revert bruhat order.\n");
+
+	for(int i = 0; i < dq->count; i++) {
+		for(doublecoset_list_t *cur = dq->cosets[i].bruhat_lower; cur; cur = cur->next) {
+			dq->lists[nlists].next = cur->to->bruhat_higher;
+			cur->to->bruhat_higher = &dq->lists[nlists];
+			cur->to->bruhat_higher->to = &dq->cosets[i];
+			nlists++;
+		}
+	}
+
+	LOG("Find opposites.\n");
+
+	for(int i = 0; i < dq->count; i++) {
+		int oppind = bitmask_reverse[~bitmask[i] & ((1<<rank) - 1)];
+		dq->cosets[i].opposite = &dq->cosets[oppind];
+	}
+
+	/*
+	printf("\n");
+	printf("digraph test123 {\n");
+
+	for(int i = 0; i < dq->count; i++) {
+		for(doublecoset_list_t *cur = dq->cosets[i].bruhat_lower; cur; cur = cur->next) {
+			printf("\"0x%02x\" -> \"0x%02x\";\n", bitmask[i], bitmask[cur->to->index]);
+		}
+	}
+
+	printf("}\n\n");
+
+	printf("Opposites:\n");
+	for(int i = 0; i < dq->count; i++)
+		printf("0x%02x <-> %d 0x%02x\n", bitmask[i], dq->cosets[i].opposite->index, bitmask[dq->cosets[i].opposite->index]);
+	printf("\n");
+	*/
+
+	long count = enumerate_balanced_thickenings(dq, 0, 0);
+	printf("Found %ld balanced ideals.\n", count);
+
+	free(bitmask);
+	free(bitmask_reverse);
+	free(dq->lists);
+	free(dq->cosets);
+	free(dq);
+
+	return 0;
+}

Makefile

@@ -4,10 +4,10 @@ HEADERS=weyl.h thickenings.h queue.h bitvec.h
 #SPECIAL_OPTIONS=-O3 -pg -funroll-loops -fno-inline
 SPECIAL_OPTIONS=-O3 -flto -funroll-loops -Winline
 
-OPTIONS=-m64 -march=native -mtune=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
+OPTIONS=-m64 -march=native -std=gnu99 -D_GNU_SOURCE $(SPECIAL_OPTIONS)
 NAME=enumerate-balanced-ideals
 
-all: enumerate graph
+all: enumerate graph D2n dominant_weights
 
 $(NAME).tar.bz2: $(NAME) $(HEADERS) enumerate.c weyl.c thickenings.c
 	tar cjhf $(NAME).tar.bz2 $(NAME)/enumerate.c $(NAME)/weyl.c $(NAME)/thickenings.c $(NAME)/weyl.h $(NAME)/thickenings.h $(NAME)/queue.h $(NAME)/bitvec.h $(NAME)/Makefile $(NAME)/graph.c
@@ -21,17 +21,29 @@ enumerate: enumerate.o weyl.o thickenings.o
 graph: graph.o weyl.o
 	gcc $(OPTIONS) -o graph graph.o weyl.o
 
+D2n: D2n.o weyl.o thickenings.o
+	gcc $(OPTIONS) -o D2n D2n.o weyl.o thickenings.o
+
+dominant_weights: dominant_weights.o weyl.o thickenings.o
+	gcc $(OPTIONS) -o dominant_weights dominant_weights.o weyl.o thickenings.o -lcdd
+
 enumerate.o: enumerate.c $(HEADERS)
 	gcc $(OPTIONS) -c enumerate.c
 
 thickenings.o: thickenings.c $(HEADERS)
 	gcc $(OPTIONS) -c thickenings.c
 
+D2n.o: D2n.c $(HEADERS)
+	gcc $(OPTIONS) -c D2n.c
+
 weyl.o: weyl.c $(HEADERS)
 	gcc $(OPTIONS) -c weyl.c
 
 graph.o: graph.c $(HEADERS)
 	gcc $(OPTIONS) -c graph.c
 
+dominant_weights.o: dominant_weights.c $(HEADERS)
+	gcc $(OPTIONS) -c dominant_weights.c
+
 clean:
-	rm -f enumerate graph thickenings.o weyl.o enumerate.o graph.o $(NAME) $(NAME).tar.bz2
+	rm -f enumerate graph D2n dominant_weights thickenings.o weyl.o enumerate.o graph.o D2n.o dominant_weights.o $(NAME) $(NAME).tar.bz2

dominant_weights.c

@@ -0,0 +1,180 @@
+#include <stdio.h>
+#include <string.h>
+#include <malloc.h>
+#include <setoper.h>
+#include <cdd.h>
+
+#include "weyl.h"
+#include "thickenings.h"
+
+doublequotient_t *dq;
+double *vector;
+char buf[1000], buf2[1000];
+
+typedef struct {
+	dd_MatrixPtr M;
+	dd_LPSolutionPtr lps;
+} info_t;
+
+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++)
+      str[j] = e->word[j] + 'a';
+    str[e->wordlength] = 0;
+  }
+
+  return str;
+}
+
+static char* print_vector(double *v, int rank, char *buf)
+{
+	int written = 0;
+
+	for(int i = 0; i < rank; i++) {
+		written += sprintf(buf+written, "%.4f", v[i]);
+		if(i != rank -1) {
+			sprintf(buf+written, ", ");
+			written += 2;
+		}
+	}
+
+	return buf;
+}
+
+static char* print_vector_ptr(dd_Arow v, int rank, char *buf)
+{
+	int written = 0;
+
+	for(int i = 0; i < rank; i++) {
+		written += sprintf(buf+written, "%.4f", v[i][0]);
+		if(i != rank -1) {
+			sprintf(buf+written, ", ");
+			written += 2;
+		}
+	}
+
+	return buf;
+}
+
+static void balanced_thickening_callback(const bitvec_t *pos, int size, const enumeration_info_t *ei)
+{
+	int bit, funcbit, sign;
+	dd_rowset ImL, Lbasis;
+	dd_ErrorType err = dd_NoError;
+	info_t *info = (info_t*)ei->callback_data;
+	int rank = weyl_rank(dq->type);
+
+	for(int i = 0; i < size; i++) {
+		bit = i < size/2 ? bv_get_bit(pos, i) : !bv_get_bit(pos, size - 1 - i);
+		sign = bit ? 1 : -1;
+		info->M->matrix[i][0][0] = 0.0;
+		for(int j = 0; j < rank; j++)
+			info->M->matrix[i][j+1][0] = sign*vector[rank*i+j];
+//		printf("0 %.2f %.2f %.2f %d %d\n", sign*vector[3*i], sign*vector[3*i+1], sign*vector[3*i+2], bit, funcbit);
+    }
+
+	for(int i = 0; i < rank; i++) {
+		info->M->matrix[i+size][0][0] = 0.0;
+		for(int j = 0; j < rank; j++) {
+			if(i == j)
+				info->M->matrix[i+size][j+1][0] = 1.0;
+			else
+				info->M->matrix[i+size][j+1][0] = 0.0;
+		}
+	}
+
+//	dd_WriteMatrix(stdout, info->M);
+
+	dd_FindRelativeInterior(info->M, &ImL, &Lbasis, &(info->lps), &err);
+
+	if(set_card(Lbasis) != 0)
+		printf("codim = %ld\n", set_card(Lbasis));
+	else
+		printf("weight = (%s)\n", print_vector_ptr(info->lps->sol, rank + 1, buf));
+
+	dd_FreeLPSolution(info->lps);
+
+//	printf("\n");
+}
+
+static int apply_reflection(double *in, double *out, int rank, int reflection)
+{
+	memcpy(out, in, rank*sizeof(double));
+	/*
+	out[reflection] *= -1;
+	if(reflection != 0)
+		out[reflection-1] += in[reflection];
+	if(reflection != rank-1)
+		out[reflection+1] += in[reflection];
+	*/
+	out[reflection] *= -1;
+	if(reflection != 0)
+		out[reflection] += in[reflection-1];
+	if(reflection != rank-1)
+		out[reflection] += in[reflection+1];
+
+}
+
+int main(int argc, char *argv[])
+{
+	semisimple_type_t type;
+	simple_type_t simple;
+	info_t info;
+
+	type.n = 1;
+	type.factors = &simple;
+	simple.series = 'A';
+	simple.rank = 3;
+
+	dq = weyl_generate_bruhat(type, 0, 0);
+
+	int order = weyl_order(type);
+	int rank = weyl_rank(type);
+	weylgroup_element_t *group = dq->group;
+	vector = malloc(weyl_order(type)*weyl_rank(type)*sizeof(double));
+
+	for(int i = 0; i < rank; i++)
+		vector[i] = 0.0;
+
+	for(int i = 0; 2*i < rank; i++) {
+		for(int j = i; j < rank - i; j++) {
+			vector[j] += rank - 2*i;
+		}
+	}
+
+	/*
+	vector[0] = 4.1;
+	vector[1] = 6.1;
+	vector[2] = 5.9;
+	vector[3] = 3.9;
+	*/
+
+	printf("regular element: (%s)\n", print_vector(vector, rank, buf));
+
+	for(int i = 0; i < order; i++) {
+		printf("%s (%s)\n", alphabetize(&group[i], buf), print_vector(vector + rank*i, rank, buf2));
+
+		for(int j = 0; j < rank; j++)
+			if(group[i].left[j]->wordlength > group[i].wordlength)
+				apply_reflection(&vector[rank*i], &vector[rank*group[i].left[j]->index], rank, j);
+	}
+
+	printf("\n");
+
+	dd_set_global_constants();
+	info.M = dd_CreateMatrix(order+rank, rank+1);
+	info.M->representation = dd_Inequality;
+	info.M->numbtype = dd_Real;
+	info.M->objective = dd_LPmax;
+
+	enumerate_balanced_thickenings(dq, balanced_thickening_callback, &info);
+
+	dd_FreeMatrix(info.M);
+	weyl_destroy_bruhat(dq);
+	free(vector);
+
+	return 0;
+}

thickenings.c

@@ -110,7 +110,7 @@ static void generate_principal_ideals(doublequotient_t *dq, bitvec_t *pos, bitve
       }
 
 #ifdef _DEBUG
-    if(is_slim[i]) {
+    if(is_slim[i] && dq->cosets[0].min) { // sometimes we don't want to define min and max
 	    fprintf(stderr, " ids: [0");
 	    for(int j = 1; j < size; j++)
 		    if(principal[j])
@@ -160,7 +160,7 @@ long enumerate_balanced_thickenings(doublequotient_t *dq, enumeration_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 < dq->count; i++)
-    if(dq->cosets[i].opposite->min->id == dq->cosets[i].min->id)
+	  if(dq->cosets[i].opposite->index == i)
       return 0;
 
   // we can only handle bitvectors up to BV_BLOCKSIZE*BV_RANK bits, but we only store half of the weyl group