colex ranking/unranking

+uses libgmp for bignum calculations. This dependency isn't very strong
as we're in most cases compiling with gcc, which itself relies on gmp.

autogen.sh

@@ -24,7 +24,7 @@ echo "noinst_HEADERS += `find lib/ -type f -name \*.h |tr \"\n\" \" \" `" >>$OUT
 echo "libcodecrypt_la_CPPFLAGS = -I\$(srcdir)/lib/ ${COMMON_CPPFLAGS}" >>$OUT
 echo "libcodecrypt_la_CFLAGS = ${COMMON_CFLAGS}" >>$OUT
 echo "libcodecrypt_la_LDFLAGS = ${COMMON_LDFLAGS}" >>$OUT
-#echo "libcodecrypt_la_LDADD = ${COMMON_LDADD} " >>$OUT
+echo "libcodecrypt_la_LIBADD = -lgmp ${COMMON_LDADD} " >>$OUT
 [ -f "lib/Makefile.am.extra" ] && cat "lib/Makefile.am.extra" >>$OUT
 
 echo "bin_PROGRAMS = ccr" >>$OUT

include/codecrypt.h

@@ -23,7 +23,7 @@ typedef unsigned int uint;
 
 /*
  * vector over GF(2). We rely on STL's vector<bool> == bit_vector
- * specialization for efficiency.
+ * specialization for space efficiency.
  */
 class polynomial;
 class gf2m;
@@ -36,8 +36,12 @@ public:
 	void add (const bvector&);
 	bool operator* (const bvector&); //dot product
 	bool zero() const;
+
 	void to_poly (polynomial&, gf2m&);
 	void from_poly (const polynomial&, gf2m&);
+
+	void colex_rank (bvector&);
+	void colex_unrank (bvector&, uint n, uint k);
 };
 
 /*

lib/bvector.cpp

@@ -47,3 +47,115 @@ void bvector::from_poly (const polynomial&r, gf2m&fld)
 	for (uint i = 0; i < size(); ++i)
 		item (i) = (r[i/fld.m] >> (i % fld.m) ) & 1;
 }
+
+/*
+ * utility colex (un)ranking for niederreiter and workalikes.
+ * see Ruskey's Combinatorial Generation, algorithm 4.10
+ *
+ * TODO use (external) cache for combination numbers to speed this up.
+ */
+
+#include <gmp.h>
+
+static void combination_number (uint n, uint k, mpz_t& r)
+{
+	mpz_t t;
+	if (k > n) {
+		mpz_set_ui (r, 0);
+		return;
+	}
+
+	if (k * 2 > n) k = n - k;
+
+	mpz_set_ui (r, 1);
+	mpz_init (t);
+
+	//upper part n*(n-1)*(n-2)*...*(n-k+1)
+	for (uint i = n; i > n - k; --i) {
+		mpz_swap (t, r);
+		mpz_mul_ui (r, t, i);
+	}
+	//lower part (div k!)
+	for (uint i = k; i > 1; --i) {
+		mpz_swap (t, r);
+		mpz_tdiv_q_ui (r, t, i);
+	}
+
+	mpz_clear (t);
+}
+
+static void bvector_to_mpz (bvector&v, mpz_t&r)
+{
+	mpz_set_ui (r, 0);
+	mpz_realloc2 (r, v.size() );
+	for (uint i = 0; i < v.size(); ++i)
+		if (v[i])
+			mpz_setbit (r, i);
+		else	mpz_clrbit (r, i);
+}
+
+static void mpz_to_bvector (mpz_t&x, bvector&r)
+{
+	r.resize (mpz_sizeinbase (x, 2) );
+	for (uint i = 0; i < r.size(); ++i)
+		r[i] = mpz_tstbit (x, i);
+}
+
+void bvector::colex_rank (bvector&r)
+{
+	mpz_t res, t, t2;
+	mpz_init_set_ui (res, 0);
+	mpz_init (t);
+	mpz_init (t2);
+
+	uint i, j;
+	j = 1;
+	for (i = 0; i < size(); ++i)
+		if (item (i) ) {
+			combination_number (i, j, t);
+			mpz_swap (t2, res);
+			mpz_add (res, t, t2);
+			++j;
+		}
+
+	mpz_to_bvector (res, r);
+
+	mpz_clear (t);
+	mpz_clear (t2);
+	mpz_clear (res);
+}
+
+#include <stdio.h>
+void bvector::colex_unrank (bvector&res, uint n, uint k)
+{
+	mpz_t r, t, t2;
+	mpz_init (r);
+	mpz_init (t);
+	mpz_init (t2);
+
+	bvector_to_mpz (*this, r);
+
+	res.clear();
+	res.resize (n, 0);
+
+	uint p;
+	for (uint i = k; i > 0; --i) {
+		p = i;
+		for (;;) {
+			combination_number (p, i, t);
+
+			if (mpz_cmp (t, r) > 0) break;
+			++p;
+		}
+
+		combination_number (p - 1, i, t);
+		mpz_swap (t2, r);
+		mpz_sub (r, t2, t);
+		if (p > n) continue; //overflow protection
+		res[p-1] = 1;
+	}
+
+	mpz_clear (r);
+	mpz_clear (t);
+	mpz_clear (t2);
+}