codecrypt/src/bvector.cpp

/*
 * This file is part of Codecrypt.
 *
 * Codecrypt is free software: you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or (at
 * your option) any later version.
 *
 * Codecrypt is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with Codecrypt. If not, see <http://www.gnu.org/licenses/>.
 */

#include "bvector.h"
#include "gf2m.h"
#include "polynomial.h"

uint bvector::hamming_weight()
{
	uint r = 0;
	for (uint i = 0; i < size(); ++i) if ( (*this) [i]) ++r;
	return r;
}

void bvector::add (const bvector&a)
{
	if (a.size() > size() ) resize (a.size(), 0);
	for (uint i = 0; i < a.size(); ++i)
		item (i) = item (i) ^ a[i];
}

void bvector::add_range (const bvector&a, uint b, uint e)
{
	if (e > size() ) resize (e, 0);
	for (uint i = b; i < e; ++i)
		item (i) = item (i) ^ a[i];
}

void bvector::add_offset (const bvector&a, uint offset)
{
	if (offset + a.size() > size() ) resize (offset + a.size(), 0);
	for (uint i = 0; i < a.size(); ++i)
		item (offset + i) = item (offset + i) ^ a[i];
}

void bvector::set_block (const bvector&a, uint offset)
{
	if (offset + a.size() > size() ) resize (offset + a.size(), 0);
	for (uint i = 0; i < a.size(); ++i)
		item (offset + i) = a[i];
}

void bvector::get_block (uint offset, uint bs, bvector&out) const
{
	if (offset + bs > size() ) return;
	out.resize (bs);
	for (uint i = 0; i < bs; ++i) out[i] = item (offset + i);
}

bool bvector::operator* (const bvector&a)
{
	bool r = 0;
	uint s = size(), i;
	if (s > a.size() ) s = a.size();
	for (i = 0; i < s; ++i) r ^= (item (i) &a[i]);
	return r;
}

bool bvector::zero() const
{
	for (uint i = 0; i < size(); ++i) if (item (i) ) return false;
	return true;
}

void bvector::to_poly (polynomial&r, gf2m&fld) const
{
	r.clear();
	if (size() % fld.m) return; //impossible
	r.resize (size() / fld.m, 0);
	for (uint i = 0; i < size(); ++i)
		if (item (i) ) r[i / fld.m] |= (1 << (i % fld.m) );
}

void bvector::from_poly (const polynomial&r, gf2m&fld)
{
	clear();
	resize (r.size() *fld.m, 0);
	for (uint i = 0; i < size(); ++i)
		item (i) = (r[i / fld.m] >> (i % fld.m) ) & 1;
}

void bvector::to_poly_cotrace (polynomial&r, gf2m&fld) const
{
	r.clear();
	if (size() % fld.m) return; //impossible
	uint s = size() / fld.m;
	r.resize (s, 0);
	for (uint i = 0; i < size(); ++i)
		if (item (i) ) r[i % s] |= (1 << (i / s) );
}

void bvector::from_poly_cotrace (const polynomial&r, gf2m&fld)
{
	clear();
	uint s = r.size();
	resize (s * fld.m, 0);
	for (uint i = 0; i < size(); ++i)
		item (i) = (r[i % s] >> (i / s) ) & 1;
}

bool bvector::to_string (std::string& out) const
{
	if (size() & 0x7) return false;

	out.clear();
	out.resize (size() >> 3, 0);

	for (uint i = 0; i < size(); ++i)
		if (item (i) ) out[i >> 3] |= (1 << (i & 0x7) );

	return true;
}

void bvector::from_string (const std::string&in)
{
	clear();
	resize (in.length() << 3);

	for (uint i = 0; i < size(); ++i)
		item (i) = (in[i >> 3] >> (i & 0x7) ) & 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 (const 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) const
{
	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);
}

void bvector::colex_unrank (bvector&res, uint n, uint k) const
{
	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);

	for (uint i = k; i > 0; --i) {
		/* Original code:

		uint p = i;
		for (;;) {
			combination_number (p, i, t);

			if (mpz_cmp (t, r) > 0) break;
			++p;
		}

		 * ...that kindof lacks speed. We're actually trying to find
		 * the smallest value of p for which comb(p,i)>r.
		 *
		 * Computing all combination numbers is KIND OF slow and cache
		 * (as suggested by Barreto and others) doesn't really make big
		 * difference here (we're usually doing only one or two runs of
		 * this stuff in one program run). Storing about 50megs of
		 * precalculated combination numbers is weird as well.
		 *
		 * Therefore, with the knowledge that i <= p <= n, we're
		 * halving the search interval as usual.
		 */

		uint p, a = i, b = n + 1;
		while (a < b) {
			p = (a + b) / 2;

			combination_number (p, i, t);
			if (mpz_cmp (t, r) > 0) b = p;
			else a = p + 1;
		}

		combination_number (b - 1, i, t);
		mpz_swap (t2, r);
		mpz_sub (r, t2, t);

		//overflow protection (result's wrong anyway now)
		if (b > n) continue;
		res[b - 1] = 1;
	}

	mpz_clear (r);
	mpz_clear (t);
	mpz_clear (t2);
}