From 1811cf143ea98ddf6a41f5dfc72d0fc9e28de48f Mon Sep 17 00:00:00 2001 From: Florian Stecker Date: Wed, 16 Nov 2016 22:59:35 +0100 Subject: [PATCH] still not working --- Makefile | 2 +- bitvec.h | 43 ++++++++++++++----- thickenings.c | 114 +++++++++++++++++++++++++++++++++++++++++--------- 3 files changed, 129 insertions(+), 30 deletions(-) diff --git a/Makefile b/Makefile index d308bf4..31883f5 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ HEADERS=coxeter.h thickenings.h queue.h bitvec.h #OPTIONS=-O3 -m64 -march=native -flto -funroll-loops -std=gnu99 -OPTIONS=-O0 -g -std=gnu99 +OPTIONS=-m64 -march=native -O0 -g -std=gnu99 #OPTIONS=-O3 -m64 -march=native -funroll-loops -fno-inline -std=gnu99 -pg all: generate process diff --git a/bitvec.h b/bitvec.h index 3d4ec6e..8df96ec 100644 --- a/bitvec.h +++ b/bitvec.h @@ -27,6 +27,10 @@ /* For _mm_popcnt_u64 */ #include +// FIRSTBITS(n) only yields useful result when 0 <= n < 64 +#define FIRSTBITS(n) ((1l << (n)) - 1l) +#define ALLBITS ((uint64_t)-1) + typedef struct { uint64_t v[BV_QWORD_RANK]; } bitvec_t; @@ -129,6 +133,8 @@ static inline void bv_print(FILE *f, const bitvec_t *x, int len) { for(int i = 0; i < len; i++) { fputc(bv_get_bit(x, i) ? '1' : '0', f); + if(i % 64 == 63) + fputc('-',f); } } @@ -161,22 +167,39 @@ static inline int bv_full(const bitvec_t *x, int len) { int i; for(i = 0; i < len/64; i++) - if(x->v[i] != (uint64_t)-1) + if(x->v[i] != ALLBITS) return 0; - return (x->v[i] & ((1<<(len%64))-1)) == ((1<<(len%64))-1); + return (x->v[i] & FIRSTBITS(len%64)) == FIRSTBITS(len%64); } -// most sigificant set bit in bits 0...start-1 -static inline int bv_last_zero_before(const bitvec_t *x, int start) +// set bits in range start...end (including start and excluding end) +static inline void bv_set_range(bitvec_t *x, int start, int end) { - // naive implementation - int i; - for(i = start - 1; i >= 0; i--) - if(!bv_get_bit(x, i)) - break; + if(start/64 == end/64) + x->v[start/64] |= ~FIRSTBITS(start % 64) & FIRSTBITS(end % 64); + else { + x->v[start/64] |= ~FIRSTBITS(start % 64); + for(int i = start/64 + 1; i < end/64; i++) + x->v[i] = ALLBITS; + x->v[end/64] |= FIRSTBITS(end % 64); + } +} - return i; +// find least significant 0 bit starting from position start (included) +static inline int bv_next_zero(const bitvec_t *x, int start) +{ + int position, i; + + position = ffsll(~(x->v[start/64] | FIRSTBITS(start % 64))); + + if(position) // found zero in same chunk? + return (start/64)%64 + position - 1; + + for(i = start/64 + 1, position = 0; i < BV_QWORD_RANK && position == 0; i++) + position = ffsll(~x->v[i]); + + return position ? (i - 1)*64 + position - 1 : BV_QWORD_RANK; } /* static inline int bv_lowest_bit_set(bitvec_t x) */ diff --git a/thickenings.c b/thickenings.c index 0b65cd1..b8ccbcd 100644 --- a/thickenings.c +++ b/thickenings.c @@ -665,20 +665,93 @@ static inline void output_thickening(const enumeration_info_t *info, signed char how can we handle negatives right to head?! */ +/* + + ok this is screwed up, let's start over: + + pos and neg are bitvectors of size info.size/2 + they stand for the first (shortest) info.size/2 elements of the weyl group + the least siginficant bit is the identity + + at the beginning they are cleared, meaning everything is undecided + */ + + +// enumerate all balanced ideals which contain everything in pos and everything strictly between known_until and next_neg, and which do not contain everything in neg, as well as next_neg +// assuming all bits up to known_until are set either in pos or in neg, but not in both +// returns number of found balanced ideals +// next_neg can be info.size/2; in that case, everything between known_until and info.size/2 is required to be in the ideal, but it does not mean that next_neg is really not contained in the ideal +// next_neg must be strictly greater than known_until, and less or equal to info.size/2 +static long enumerate_tree(const enumeration_info_t *info, const bitvec_t *pos, const bitvec_t *neg, int first_unknown, int next_neg) +{ + bitvec_t newpos, newneg, known; + int next_next_neg; + long count = 0; + + // the omission of next_neg means inclusion of info->size - 1 - next_neg + // add its principal ideal to pos and the opposite to neg + if(next_neg != info->size/2) { + bv_union(&info->principal_pos[info->size - 1 - next_neg], pos, &newpos); + bv_union(&info->principal_neg[info->size - 1 - next_neg], neg, &newneg); + } + + // add the range from first_unknown to next_neg to newpos + bv_set_range(&newpos, first_unknown, next_neg); // including the start, excluding end + + bitvec_t test1, test2; + bv_print(stderr, &newpos, info->size/2); + fprintf(stderr, " "); + bv_print(stderr, &newneg, info->size/2); + fprintf(stderr, " "); + bv_union(&newpos, &newneg, &test1); + bv_print(stderr, &test1, info->size/2); + fprintf(stderr, " %d\n", bv_disjoint(&newpos, &newneg)); + + + // check if this leads to any conflicts (equivalently, violates slimness) + if(!bv_disjoint(&newpos, &newneg)) + return 0; + + // what do we know so far? + bv_union(&newpos, &newneg, &known); + + // do we know everything already? we have a balanced ideal then + if(bv_full(&known, info->size/2)) { + fprintf(stderr, "Found balanced ideal: "); + bv_print(stderr, &newpos, info->size/2); + fprintf(stderr, " "); + bv_print(stderr, &newneg, info->size/2); + fprintf(stderr, "\n"); + + return 1; + } + + next_next_neg = next_neg; + while(next_next_neg < info->size/2) { + next_next_neg = bv_next_zero(&known, next_next_neg + 1); // this could return info->size/2, but that's fine for enumerate_tree + count += enumerate_tree(info, &newpos, &newneg, next_neg + 1, next_next_neg); + } +} + +/* static long enumerate_tree(const enumeration_info_t *info, const bitvec_t *pos, const bitvec_t *neg, int head, int oldhead) { bitvec_t newpos, newneg, known; int newhead; long count = 0; - bv_union(&info->principal_pos[head], pos, &newpos); - bv_union(&info->principal_neg[head], neg, &newneg); - - for(int i = oldhead - 1; i > head; i--) { - bv_union(&info->principal_pos[info->size-1-i], &newpos, &newpos); - bv_union(&info->principal_neg[info->size-1-i], &newneg, &newneg); + if(head < info->size/2) { + fprintf(stderr, "head is %d!\n", head); } + if(head != -1) { + bv_union(&info->principal_pos[head], pos, &newpos); + bv_union(&info->principal_neg[head], neg, &newneg); + } + + for(int i = oldhead - 1; i > head && i >= info->size/2; i--) + bv_set_bit(&newpos, info->size - 1 - i); + if(!bv_disjoint(&newpos, &newneg)) return 0; @@ -696,16 +769,16 @@ static long enumerate_tree(const enumeration_info_t *info, const bitvec_t *pos, return 1; } - newhead = bv_last_zero_before(&known, head); - while(newhead >= 0) { + newhead = head; + while(newhead >= info->size/2) { + newhead = bv_last_zero_before(&known, newhead); // this will return -1 at some point, which stands for the possibility that all remaining choices are "no" count += enumerate_tree(info, &newpos, &newneg, newhead, head); - newhead = bv_last_zero_before(&known, newhead); } return count; } -/* + static long enumerate_tree(const enumeration_info_t *info, signed char *level, int current_level, int head) { @@ -822,19 +895,22 @@ long enumerate_balanced_thickenings(semisimple_type_t type, node_t *graph, int s fprintf(stderr, "\n"); } + // test these functions a little + bitvec_t test; + bv_clear(&test); + bv_set_range(&test, 0, 4); + bv_print(stderr, &test, 200); + fprintf(stderr, " %d\n", bv_full(&test, 6)); + + // don't do enumeration + // return count; + // enumerate balanced ideals bitvec_t pos, neg; bv_clear(&pos); bv_clear(&neg); - for(int i = info.size - 1; i >= 0; i--) - count += enumerate_tree(&info, &pos, &neg, i, info.size); - - /* old algorithm: - level = (signed char*)malloc(info.size*sizeof(int)); - memset(level, 0, info.size*sizeof(int)); - - for(int i = info.size - 1; i >= 0; i--) - count += enumerate_tree(&info, level, 1, i); */ + for(int i = 0; i <= info.size/2; i++) + count += enumerate_tree(&info, &pos, &neg, 0, i); free(principal_pos); free(principal_neg);