Rename library "chilbert"

10 files changed, 2286 insertions, 0 deletions

diff --git a/chilbert/Algorithm.hpp b/chilbert/Algorithm.hpp
new file mode 100644
index 0000000..81bd050
--- /dev/null
+++ b/chilbert/Algorithm.hpp
@@ -0,0 +1,476 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _ALGORITHM_HPP_
+#define _ALGORITHM_HPP_
+
+
+#include "chilbert/BigBitVec.hpp"
+#include "chilbert/GetLocation.hpp"
+#include "chilbert/SetLocation.hpp"
+#include "chilbert/GetBits.hpp"
+#include "chilbert/SetBits.hpp"
+#include "chilbert/GrayCodeRank.hpp"
+#include <string.h>
+
+
+// Templated Hilbert functions.
+// P - is the class used to represent each dimension
+//     of a multidimensional point.
+// H - is the class used to represent the point as a Hilbert
+//     index.
+// I - is the class used to represent interim variables.
+//     Needs to have as many bits as there are dimensions.
+//
+// In general, each of P,H,I should be a FixBitVec if they
+// fit in that fixed precision.  Otherwise, use a BigBitVec.
+// Whatever you use, they must all be of consistent underlying
+// storage types.
+
+
+// The dimension across which the Hilbert curve travels
+// principally.
+// D0 is d + 1, where d is the actual dimension you want to
+// walk across.
+// MUST HAVE 0 <= D0 < n
+#define D0	1
+
+namespace chilbert
+{
+	// 'Transforms' a point.
+	template<class I>
+	inline
+	void
+	transform(
+		const I &e,
+		int d,
+		int n,
+		I &a
+		)
+	{
+		a ^= e;
+		a.rotr( d, n );//#D d+1, n );
+		return;
+	}
+
+	// Inverse 'transforms' a point.
+	template<class I>
+	inline
+	void
+	transformInv(
+		const I &e,
+		int d,
+		int n,
+		I &a
+		)
+	{
+		a.rotl( d, n );//#D d+1, n );
+		a ^= e;
+		return;
+	}
+
+	// Update for method 1 (GrayCodeInv in the loop)
+	template<class I>
+	inline
+	void
+	update1(
+		const I &l,
+		const I &t,
+		const I &w,
+		int n,
+		I &e,
+		int &d
+		)
+	{
+		assert( 0 <= d && d < n );
+		e = l;
+		e.toggle( d ); //#D d == n-1 ? 0 : d+1 );
+	
+		// Update direction
+		d += 1 + t.fsb();
+		if ( d >= n ) d -= n;
+		if ( d >= n ) d -= n;
+		assert( 0 <= d && d < n );
+
+		if ( ! (w.rack() & 1) )
+			e.toggle( d == 0 ? n-1 : d-1 ); //#D d );
+
+		return;
+	}
+
+	// Update for method 2 (GrayCodeInv out of loop)
+	template<class I>
+	inline
+	void
+	update2(
+		const I &l,
+		const I &t,
+		const I &w,
+		int n,
+		I &e,
+		int &d
+		)
+	{
+		assert( 0 <= d && d < n );
+		e = l;
+		e.toggle( d );//#D d == n-1 ? 0 : d+1 );
+	
+		// Update direction
+		d += 1 + t.fsb();
+		if ( d >= n ) d -= n;
+		if ( d >= n ) d -= n;
+		assert( 0 <= d && d < n );
+
+		return;
+	}
+
+	template <class P,class H,class I>
+	inline
+	void
+	_coordsToIndex(
+		const P *p,
+		int m,
+		int n,
+		H &h,
+		int *ds = nullptr // #HACK
+		)
+	{
+		I e(n), l(n), t(n), w(n);
+		int d, i;
+		int ho = m*n;
+
+		// Initialize
+		e.reset();
+		d = D0;
+		l.reset();
+		h.reset();
+
+		// Work from MSB to LSB
+		for ( i = m-1; i >= 0; i-- )
+		{
+			// #HACK
+			if ( ds ) ds[i] = d;
+			
+			// Get corner of sub-hypercube where point lies.
+			getLocation<P,I>(p,n,i,l);
+
+			// Mirror and reflect the location.
+			// t = T_{(e,d)}(l)
+			t = l;
+			transform<I>(e,d,n,t);
+
+			w = t;
+			if ( i < m-1 )
+				w.toggle(n - 1);
+
+			// Concatenate to the index.
+			ho -= n;
+			setBits<H,I>(h,n,ho,w);
+
+			// Update the entry point and direction.
+			update2<I>(l,t,w,n,e,d);
+		}
+
+		h.grayCodeInv();
+
+		return;
+	}
+
+	// This is wrapper to the basic Hilbert curve index
+	// calculation function.  It will support fixed or
+	// arbitrary precision, templated.  Depending on the
+	// number of dimensions, it will use the most efficient
+	// representation for interim variables.
+	// Assumes h is big enough for the output (n*m bits!)
+	template<class P,class H>
+	inline
+	void
+	coordsToIndex(
+		const P *p,	// [in ] point
+		int m,		// [in ] precision of each dimension in bits
+		int n,		// [in ] number of dimensions
+		H &h		// [out] Hilbert index
+		)
+	{
+		// Intermediate variables will fit in fixed width?
+		if ( n <= FBV_BITS )
+			_coordsToIndex<P,H,CFixBitVec>(p,m,n,h);
+		// Otherwise, they must be BigBitVecs.
+		else
+			_coordsToIndex<P,H,CBigBitVec>(p,m,n,h);
+
+		return;
+	}
+
+
+	template <class P,class H,class I>
+	inline
+	void
+	_indexToCoords(
+		P *p,
+		int m,
+		int n,
+		const H &h
+		)
+	{
+		I e(n), l(n), t(n), w(n);
+		int d, i, j, ho;
+
+		// Initialize
+		e.reset();
+		d = D0;
+		l.reset();
+		for ( j = 0; j < n; j++ )
+			p[j] = 0;
+
+		ho = m*n;
+		
+		// Work from MSB to LSB
+		for ( i = m-1; i >= 0; i-- )
+		{
+			// Get the Hilbert index bits
+			ho -= n;
+			getBits<H,I>(h,n,ho,w);
+
+			// t = GrayCode(w)
+			t = w;
+			t.grayCode();
+
+			// Reverse the transform
+			// l = T^{-1}_{(e,d)}(t)
+			l = t;
+			transformInv<I>(e,d,n,l);
+
+			// Distribute these bits
+			// to the coordinates.
+			setLocation<P,I>(p,n,i,l);
+
+			// Update the entry point and direction.
+			update1<I>(l,t,w,n,e,d);
+		}
+
+		return;
+	}
+
+	// This is wrapper to the basic Hilbert curve inverse
+	// index function.  It will support fixed or
+	// arbitrary precision, templated.  Depending on the
+	// number of dimensions, it will use the most efficient
+	// representation for interim variables.
+	// Assumes each entry of p is big enough to hold the
+	// appropriate variable.
+	template<class P,class H>
+	inline
+	void
+	indexToCoords(
+		P *p,				// [out] point
+		int m,			// [in ] precision of each dimension in bits
+		int n,			// [in ] number of dimensions
+		const H &h	// [out] Hilbert index
+		)
+	{
+		// Intermediate variables will fit in fixed width?
+		if ( n <= FBV_BITS )
+			_indexToCoords<P,H,CFixBitVec>(p,m,n,h);
+		// Otherwise, they must be BigBitVecs.
+		else
+			_indexToCoords<P,H,CBigBitVec>(p,m,n,h);
+
+		return;
+	}
+
+	template <class P,class HC,class I>
+	inline
+	void
+	_coordsToCompactIndex(
+		const P *p,
+		const int *ms,
+		int n,
+		HC &hc,
+		int M = 0,
+		int m = 0
+		)
+	{
+		int i, mn;
+		int *ds;
+		
+		// Get total precision and max precision
+		// if not supplied
+		if ( M == 0 || m == 0 )
+		{
+			M = m = 0;
+			for ( i = 0; i < n; i++ )
+			{
+				if ( ms[i] > m ) m = ms[i];
+				M += ms[i];
+			}
+		}
+
+		mn = m*n;
+
+		// If we could avoid allocation altogether (ie: have a
+		// fixed buffer allocated on the stack) then this increases
+		// speed by a bit (4% when n=4, m=20)
+		ds = new int [ m ];
+
+		if ( mn > FBV_BITS )
+		{
+			CBigBitVec h(mn);
+			_coordsToIndex<P,CBigBitVec,I>(p,m,n,h,ds);
+			compactIndex<CBigBitVec,HC>(ms,ds,n,m,h,hc);
+		}
+		else
+		{
+			CFixBitVec h;
+			_coordsToIndex<P,CFixBitVec,I>(p,m,n,h,ds);
+			compactIndex<CFixBitVec,HC>(ms,ds,n,m,h,hc);
+		}
+
+		delete [] ds;
+
+		return;
+	}
+
+	// This is wrapper to the basic Hilbert curve index
+	// calculation function.  It will support fixed or
+	// arbitrary precision, templated.  Depending on the
+	// number of dimensions, it will use the most efficient
+	// representation for interim variables.
+	// Assumes h is big enough for the output (n*m bits!)
+	template<class P,class HC>
+	inline
+	void
+	coordsToCompactIndex(
+		const P *p,		// [in ] point
+		const int *ms,// [in ] precision of each dimension in bits
+		int n,				// [in ] number of dimensions
+		HC &hc,					// [out] Hilbert index
+		int M = 0,
+		int m = 0
+		)
+	{
+		// Intermediate variables will fit in fixed width?
+		if ( n <= FBV_BITS )
+			_coordsToCompactIndex<P,HC,CFixBitVec>(p,ms,n,hc,M,m);
+		// Otherwise, they must be BigBitVecs.
+		else
+			_coordsToCompactIndex<P,HC,CBigBitVec>(p,ms,n,hc,M,m);
+
+		return;
+	}
+
+	template <class P,class HC,class I>
+	inline
+	void
+	_compactIndexToCoords(
+		P *p,
+		const int *ms,
+		int n,
+		const HC &hc,
+		int M = 0,
+		int m = 0
+		)
+	{
+		I e(n), l(n), t(n), w(n), r(n), mask(n), ptrn(n);
+		int d, i, j, b;
+
+		// Get total precision and max precision
+		// if not supplied
+		if ( M == 0 || m == 0 )
+		{
+			M = m = 0;
+			for ( i = 0; i < n; i++ )
+			{
+				if ( ms[i] > m ) m = ms[i];
+				M += ms[i];
+			}
+		}
+		
+		// Initialize
+		e.reset();
+		d = D0;
+		l.reset();
+		for ( j = 0; j < n; j++ )
+			p[j].reset();
+		
+		// Work from MSB to LSB
+		for ( i = m-1; i >= 0; i-- )
+		{
+			// Get the mask and ptrn
+			extractMask<I>(ms,n,d,i,mask,b);
+			ptrn = e;
+			ptrn.rotr(d,n);//#D ptrn.Rotr(d+1,n);
+
+			// Get the Hilbert index bits
+			M -= b;
+			r.reset(); // GetBits doesn't do this
+			getBits<HC,I>(hc,b,M,r);
+
+			// w = GrayCodeRankInv(r)
+			// t = GrayCode(w)
+			grayCodeRankInv<I>(mask,ptrn,r,n,b,t,w);
+
+			// Reverse the transform
+			// l = T^{-1}_{(e,d)}(t)
+			l = t;
+			transformInv<I>(e,d,n,l);
+
+			// Distribute these bits
+			// to the coordinates.
+			setLocation<P,I>(p,n,i,l);
+
+			// Update the entry point and direction.
+			update1<I>(l,t,w,n,e,d);
+		}
+
+		return;
+	}
+
+	// This is wrapper to the basic Hilbert curve inverse
+	// index function.  It will support fixed or
+	// arbitrary precision, templated.  Depending on the
+	// number of dimensions, it will use the most efficient
+	// representation for interim variables.
+	// Assumes each entry of p is big enough to hold the
+	// appropriate variable.
+	template<class P,class HC>
+	inline
+	void
+	compactIndexToCoords(
+		P *p,					// [out] point
+		const int *ms,// [in ] precision of each dimension in bits
+		int n,				// [in ] number of dimensions
+		const HC &hc,		// [out] Hilbert index
+		int M = 0,
+		int m = 0
+		)
+	{
+		// Intermediate variables will fit in fixed width?
+		if ( n <= FBV_BITS )
+			_compactIndexToCoords<P,HC,CFixBitVec>(p,ms,n,hc,M,m);
+		// Otherwise, they must be BigBitVecs.
+		else
+			_compactIndexToCoords<P,HC,CBigBitVec>(p,ms,n,hc,M,m);
+
+		return;
+	}
+
+};
+
+
+#endif
diff --git a/chilbert/BigBitVec.hpp b/chilbert/BigBitVec.hpp
new file mode 100644
index 0000000..2f924b9
--- /dev/null
+++ b/chilbert/BigBitVec.hpp
@@ -0,0 +1,799 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _BIGBITVEC_HPP_
+#define _BIGBITVEC_HPP_
+
+
+#include "chilbert/FixBitVec.hpp"
+
+#include <cstring>
+
+#define BBV_MIN(a,b)		((a)<(b)?(a):(b))
+#define BBV_MAX(a,b)		((a)>(b)?(a):(b))
+#define FBVS_NEEDED(b)		((BBV_MAX(b,1)+FBV_BITS-1)/FBV_BITS)
+#define BBV_MODSPLIT(r,b,k) { b=(k); r=b/FBV_BITS; b-=r*FBV_BITS; }
+
+namespace chilbert {
+
+class CBigBitVec
+{
+public:
+
+	// Constructor, with optional number of bits.
+	CBigBitVec(
+		int iBits = FBV_BITS
+	           )
+	{
+		// Determine number of racks required.
+		m_iRacks = iBits == 0 ? 0 : FBVS_NEEDED(iBits);
+
+		// Allocate the memory.
+		if (m_iRacks > 0) {
+			m_pcRacks = new CFixBitVec[m_iRacks];
+		} else {
+			m_pcRacks = nullptr;
+		}
+			
+		return;
+	}
+		
+		
+	// Copy construct.	Creates duplicate.
+	CBigBitVec(
+		const CBigBitVec &cBBV
+	           )
+	{
+		m_iRacks = cBBV.m_iRacks;
+		m_pcRacks = new CFixBitVec[m_iRacks];
+
+		// Copy the rack values.
+		/*int i;
+		  for ( i = 0; i < m_iRacks; i++ )
+		  m_pcRacks[i] = cBBV.m_pcRacks[i];*/
+		if (cBBV.m_pcRacks) {
+			memcpy( static_cast<void*>(m_pcRacks),
+			        static_cast<const void*>(cBBV.m_pcRacks),
+			        sizeof(CFixBitVec)*m_iRacks );
+		}
+
+		return;
+	}
+
+	// Move construct.
+	CBigBitVec(
+		CBigBitVec &&cBBV
+	           )
+	{
+		m_iRacks = cBBV.m_iRacks;
+		m_pcRacks = cBBV.m_pcRacks;
+		cBBV.m_iRacks = 0;
+		cBBV.m_pcRacks = nullptr;
+	}
+
+	// Copy constructor.
+	CBigBitVec(
+		const CFixBitVec &cFBV
+	           )
+	{
+		m_iRacks = 1;
+		m_pcRacks = new CFixBitVec[m_iRacks];
+
+		m_pcRacks[0] = cFBV;
+
+		return;
+	}
+
+
+	// Destructor
+	~CBigBitVec()
+	{
+		delete [] m_pcRacks;
+
+		return;
+	}
+
+
+	// Returns the current size in bits.
+	int
+	size() const
+	{
+		return m_iRacks*FBV_BITS;
+	}
+
+
+	// zeros the bit-vector.
+	CBigBitVec &
+	reset()
+	{
+	
+		/*int i;
+		  for ( i = 0; i < m_iRacks; i++ )
+		  m_pcRacks[i].reset();*/
+		memset( static_cast<void*>(m_pcRacks), 0,
+		        sizeof(CFixBitVec)*m_iRacks );
+
+		return (*this);
+	}
+
+
+	// truncates the bit-vector to a given precision in
+	// bits (zeroes MSBs without shrinking the vector)
+	CBigBitVec &
+	truncate(
+		int iBits
+	         )
+	{
+		assert( iBits >= 0 && iBits <= size() );
+		int r, b, i;
+
+		BBV_MODSPLIT(r,b,iBits);
+	
+		if ( r >= m_iRacks )
+			return (*this);
+
+		m_pcRacks[r].truncate(b);
+	
+		for ( i = r+1; i < m_iRacks; i++ )
+			m_pcRacks[i].reset();
+
+		return (*this);
+	}
+
+
+	// Assignment operator.	No resizing.
+	CBigBitVec &
+	operator=(
+		const CBigBitVec &cBBV
+	          )
+	{
+		if ( m_iRacks < cBBV.m_iRacks )
+		{
+			/*int i;
+			  for ( i = 0; i < m_iRacks; i++ )
+			  m_pcRacks[i] = cBBV.m_pcRacks[i];*/
+			memcpy( static_cast<void*>(m_pcRacks),
+			        static_cast<const void*>(cBBV.m_pcRacks),
+			        sizeof(CFixBitVec)*m_iRacks );
+		}
+		else
+		{
+			/*int i;
+			  for ( i = 0; i < cBBV.m_iRacks; i++ )
+			  m_pcRacks[i] = cBBV.m_pcRacks[i];
+			  for ( ; i < m_iRacks; i++ )
+			  m_pcRacks[i].reset();*/
+			if (m_pcRacks) {
+				if (cBBV.m_pcRacks) {
+				memcpy( static_cast<void*>(m_pcRacks),
+				        static_cast<const void*>(cBBV.m_pcRacks),
+				        sizeof(CFixBitVec)*cBBV.m_iRacks );
+				}
+				memset( static_cast<void*>(m_pcRacks+cBBV.m_iRacks),
+				        0, sizeof(CFixBitVec)*(m_iRacks-cBBV.m_iRacks) );
+			}
+		}
+
+		return (*this);
+	}
+	CBigBitVec &
+	operator=(
+		CBigBitVec &&cBBV
+	          )
+	{
+		if (&cBBV != this) {
+			m_iRacks = cBBV.m_iRacks;
+			m_pcRacks = cBBV.m_pcRacks;
+			cBBV.m_iRacks = 0;
+			cBBV.m_pcRacks = nullptr;
+		}
+
+		return (*this);
+	}
+	CBigBitVec &
+	operator=(
+		const CFixBitVec &cFBV
+	          )
+	{
+		m_pcRacks[0] = cFBV;
+		/*int i;
+		  for ( i = 1; i < m_iRacks; i++ )
+		  m_pcRacks[i].reset();*/
+		memset( static_cast<void*>(m_pcRacks+1),
+		        0, sizeof(CFixBitVec)*(m_iRacks-1) );
+		return (*this);
+	}
+	CBigBitVec &
+	operator=(
+		FBV_UINT j
+	          )
+	{
+		m_pcRacks[0] = j;
+		/*int i;
+		  for ( i = 1; i < m_iRacks; i++ )
+		  m_pcRacks[i].reset();*/
+		memset( static_cast<void*>(m_pcRacks+1),
+		        0, sizeof(CFixBitVec)*(m_iRacks-1) );
+		return (*this);
+	}
+
+	// Returns the value of the nth bit.
+	bool
+	test(
+		int iIndex
+	       ) const
+	{
+		assert( iIndex >= 0 && iIndex < size() );
+		int r, b;
+		BBV_MODSPLIT(r,b,iIndex);
+		return m_pcRacks[r].test(b);
+	}
+
+
+	// Sets the value of the nth bit.
+	CBigBitVec &
+	set(
+		int iIndex,
+		bool bBit = true
+	       )
+	{
+		assert( iIndex >= 0 && iIndex < size() );
+		int r, b;
+		BBV_MODSPLIT(r,b,iIndex);
+		m_pcRacks[r].set(b,bBit);
+		return (*this);
+	}
+
+
+	// Toggles the value of the nth bit.
+	CBigBitVec &
+	toggle(
+		int iIndex
+	          )
+	{
+		assert( iIndex >= 0 && iIndex < size() );
+		int r, b;
+		BBV_MODSPLIT(r,b,iIndex);
+		m_pcRacks[r].toggle(b);
+		return (*this);
+	}
+
+
+	// In place AND.
+	CBigBitVec &
+	operator&=(
+		const CBigBitVec &cBBV
+	           )
+	{
+		int i;
+	
+		for ( i = 0; i < BBV_MIN(m_iRacks,cBBV.m_iRacks); i++ )
+			m_pcRacks[i] &= cBBV.m_pcRacks[i];
+	
+		return (*this);
+	}
+	CBigBitVec &
+	operator&=(
+		const CFixBitVec &r
+	           )
+	{
+		m_pcRacks[0] &= r;
+		return (*this);
+	}
+	CBigBitVec &
+	operator&=(
+		FBV_UINT i
+	           )
+	{
+		m_pcRacks[0] &= i;
+		return (*this);
+	}
+
+
+	// AND operator.
+	CBigBitVec
+	operator&(
+		const CBigBitVec &cBBV
+	          ) const
+	{
+		CBigBitVec t( *this );
+		t &= cBBV;
+
+		return t;
+	}
+	CBigBitVec
+	operator&(
+		const CFixBitVec &r
+	          )
+	{
+		CBigBitVec t( *this );
+		t &= r;
+	
+		return t;
+	}
+	CBigBitVec
+	operator&(
+		FBV_UINT i
+	          )
+	{
+		CBigBitVec t( *this );
+		t &= i;
+	
+		return t;
+	}
+
+
+	// In place OR.
+	CBigBitVec &
+	operator|=(
+		const CBigBitVec &cBBV
+	           )
+	{
+		int i;
+	
+		for ( i = 0; i < BBV_MIN(m_iRacks,cBBV.m_iRacks); i++ )
+			m_pcRacks[i] |= cBBV.m_pcRacks[i];
+	
+		return (*this);
+	}
+	CBigBitVec &
+	operator|=(
+		const CFixBitVec &r
+	           )
+	{
+		m_pcRacks[0] |= r;
+		return (*this);
+	}
+	CBigBitVec &
+	operator|=(
+		FBV_UINT i
+	           )
+	{
+		m_pcRacks[0] |= i;
+		return (*this);
+	}
+
+
+	// OR operator.
+	CBigBitVec
+	operator|(
+		const CBigBitVec &cBBV
+	          ) const
+	{
+		CBigBitVec t( *this );
+		t |= cBBV;
+
+		return t;
+	}
+	CBigBitVec
+	operator|(
+		const CFixBitVec &r
+	          )
+	{
+		CBigBitVec t( *this );
+		t |= r;
+	
+		return t;
+	}
+	CBigBitVec
+	operator|(
+		FBV_UINT i
+	          )
+	{
+		CBigBitVec t( *this );
+		t |= i;
+	
+		return t;
+	}
+
+
+	// In place XOR.
+	CBigBitVec &
+	operator^=(
+		const CBigBitVec &cBBV
+	           )
+	{
+		int i;
+	
+		for ( i = 0; i < BBV_MIN(m_iRacks,cBBV.m_iRacks); i++ )
+			m_pcRacks[i] ^= cBBV.m_pcRacks[i];
+	
+		return (*this);
+	}
+	CBigBitVec &
+	operator^=(
+		const CFixBitVec &r
+	           )
+	{
+		m_pcRacks[0] ^= r;
+		return (*this);
+	}
+	CBigBitVec &
+	operator^=(
+		FBV_UINT i
+	           )
+	{
+		m_pcRacks[0] ^= i;
+		return (*this);
+	}
+
+
+	// XOR operator.
+	CBigBitVec
+	operator^(
+		const CBigBitVec &cBBV
+	          ) const
+	{
+		CBigBitVec t( *this );
+		t ^= cBBV;
+
+		return t;
+	}
+	CBigBitVec	
+	operator^(
+		const CFixBitVec &r
+	          )
+	{
+		CBigBitVec t( *this );
+		t ^= r;
+	
+		return t;
+	}
+	CBigBitVec
+	operator^(
+		FBV_UINT i
+	          )
+	{
+		CBigBitVec t( *this );
+		t ^= i;
+	
+		return t;
+	}
+
+
+	// Shift left operation, in place.
+	CBigBitVec &
+	operator<<=(
+		int iBits
+	            )
+	{
+		assert( iBits >= 0 );
+		int r, b, i;
+
+		// No shift?
+		if ( iBits == 0 )
+			return (*this);
+
+		BBV_MODSPLIT(r,b,iBits);
+
+		// All racks?
+		if ( r >= m_iRacks )
+		{
+			reset();
+			return (*this);
+		}
+
+		// Do rack shifts.
+		if ( r > 0 )
+		{
+			for ( i = m_iRacks-1; i >= r; i-- )
+				m_pcRacks[i] = m_pcRacks[i-r];
+			for ( ; i >= 0; i-- )
+				m_pcRacks[i].reset();
+		}
+
+		// Do bit shifts.
+		if ( b > 0 )
+		{
+			int bi = FBV_BITS - b;
+			for ( i = m_iRacks-1; i >= r+1; i-- )
+			{
+				m_pcRacks[i] <<= b;
+				m_pcRacks[i] |= m_pcRacks[i-1] >> bi;
+			}
+			m_pcRacks[i] <<= b;
+		}
+
+		return (*this);
+	}
+
+
+	// Shift left operation.
+	CBigBitVec
+	operator<<(
+		int iBits
+	           ) const
+	{
+		CBigBitVec t( *this );
+		t <<= iBits;
+		return t;
+	}
+
+
+	// Shift right operation, in place.
+	CBigBitVec &
+	operator>>=(
+		int iBits
+	            )
+	{
+		assert( iBits >= 0 );
+		int r, b, i;
+
+		// No shift?
+		if ( iBits == 0 )
+			return (*this);
+
+		BBV_MODSPLIT(r,b,iBits);
+
+		// All racks?
+		if ( r >= m_iRacks )
+		{
+			reset();
+			return (*this);
+		}
+
+		// Do rack shifts.
+		if ( r > 0 )
+		{
+			for ( i = 0; i < m_iRacks-r; i++ )
+				m_pcRacks[i] = m_pcRacks[i+r];
+			for ( ; i < m_iRacks; i++ )
+				m_pcRacks[i].reset();
+		}
+
+		// Do bit shifts.
+		if ( b > 0 )
+		{
+			int bi = FBV_BITS - b;
+			for ( i = 0; i < m_iRacks-r-1; i++ )
+			{
+				m_pcRacks[i] >>= b;
+				m_pcRacks[i] |= m_pcRacks[i+1] << bi;
+			}
+			m_pcRacks[i] >>= b;
+		}
+
+		return (*this);
+	}
+
+
+	// Shift right operation.
+	CBigBitVec
+	operator>>(
+		int iBits
+	           ) const
+	{
+		CBigBitVec t( *this );
+		t >>= iBits;
+		return t;
+	}
+
+
+	// Right rotation, in place.
+	CBigBitVec &
+	rotr(
+		int iBits,
+		int iWidth
+	     )
+	{
+		assert( iBits >= 0 );
+
+		// Fill in the width, if necessary.
+		if ( iWidth <= 0 )
+			iWidth = size();
+
+		// Modulo the number of bits.
+		assert( iBits < iWidth );
+		if ( iBits == 0 ) return (*this);
+
+		// Ensure we are truncated appropriately.
+		truncate(iWidth);
+
+		CBigBitVec t1( *this );
+		(*this) >>= iBits;
+		t1 <<= (iWidth-iBits);
+		(*this) |= t1;
+	
+		truncate(iWidth);
+
+		return (*this);
+	}
+
+
+	// Left rotation, in place.
+	CBigBitVec &
+	rotl(
+		int iBits,
+		int iWidth
+	     )
+	{
+		assert( iBits >= 0 );
+
+		// Fill in the width, if necessary.
+		if ( iWidth <= 0 )
+			iWidth = size();
+
+		// Modulo the number of bits.
+		assert( iBits < iWidth );
+		if ( iBits == 0 ) return (*this);
+
+		// Ensure we are truncated appropriately.
+		truncate(iWidth);
+
+		CBigBitVec t1( *this );
+		(*this) <<= iBits;
+		t1 >>= (iWidth-iBits);
+		(*this) |= t1;
+
+		truncate(iWidth);
+
+		return (*this);
+	}
+
+
+	// Returns true if the rack is zero valued.
+	bool
+	none() const
+	{
+		int i;
+		for ( i = 0; i < m_iRacks; i++ )
+			if ( !m_pcRacks[i].none() ) return false;
+		return true;
+	}
+
+
+	// Returns the index of the first set bit.
+	// (numbered 1 to n, with 0 meaning no bits were set)
+	int
+	fsb() const
+	{
+		for (int i = 0; i < m_iRacks; ++i )
+		{
+			const int j = m_pcRacks[i].fsb();
+			if ( j ) {
+				return (i * FBV_BITS) + j;
+			}
+		}
+		return 0;
+	}
+
+
+	// Gray Code
+	CBigBitVec &
+	grayCode()
+	{
+		int i;
+		FBV_UINT s = 0;
+
+		for ( i = m_iRacks-1; i >= 0; i-- )
+		{
+			FBV_UINT t = m_pcRacks[i].rack() & 1;
+			m_pcRacks[i].grayCode();
+			m_pcRacks[i] ^= (s<<(FBV_BITS-1));
+			s = t;
+		}
+
+		return (*this);
+	}
+
+
+	// Gray Code Inverse
+	CBigBitVec &
+	grayCodeInv()
+	{
+		int i;
+		FBV_UINT s = 0;
+
+		for ( i = m_iRacks-1; i >= 0; i-- )
+		{
+			m_pcRacks[i].grayCodeInv();
+			if ( s ) m_pcRacks[i].flip();
+			s = m_pcRacks[i].rack() & 1;
+		}
+		return (*this);
+	}
+
+
+	// Complement
+	CBigBitVec &
+	flip()
+	{
+		int i;
+		for ( i = 0; i < m_iRacks; i++ )
+			m_pcRacks[i].flip();
+		return (*this);
+	}
+
+
+	// Returns the first rack.
+	FBV_UINT &
+	rack()
+	{
+		return m_pcRacks[0].rack();
+	}
+	FBV_UINT
+	rack() const
+	{
+		return m_pcRacks[0].rack();
+	}
+
+
+	// Returns the racks.
+	FBV_UINT *
+	racks()
+	{
+		return reinterpret_cast<FBV_UINT*>(m_pcRacks);
+	}
+	const FBV_UINT *
+	racks() const
+	{
+		return reinterpret_cast<FBV_UINT*>(m_pcRacks);
+	}
+
+
+	// Returns the number of racks.
+	int
+	rackCount() const
+	{
+		return m_iRacks;
+	}
+
+	
+private:	
+
+	// Right rotates entire racks (in place).
+	void
+	rackRotr(
+		int k
+	         )
+	{
+		assert( 0 <= k && k < m_iRacks );
+	
+		int c, v;
+		CFixBitVec tmp;
+
+		if (k == 0) return;
+
+		c = 0;
+		for (v = 0; c < m_iRacks; v++)
+		{
+			int t = v, tp = v + k;
+			tmp = m_pcRacks[v];
+			c++;
+			while (tp != v)
+			{
+				m_pcRacks[t] = m_pcRacks[tp];
+				t = tp;
+				tp += k;
+				if (tp >= m_iRacks) tp -= m_iRacks;
+				c++;
+			}
+			m_pcRacks[t] = tmp;
+		}
+
+		return;
+	}
+	
+
+	CFixBitVec *m_pcRacks;
+	int m_iRacks;
+};
+
+} // namespace chilbert
+
+#endif
diff --git a/chilbert/FixBitVec.hpp b/chilbert/FixBitVec.hpp
new file mode 100644
index 0000000..ff31f1b
--- /dev/null
+++ b/chilbert/FixBitVec.hpp
@@ -0,0 +1,421 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _FIXBITVEC_HPP_
+#define _FIXBITVEC_HPP_
+
+#include "chilbert/Operations.hpp"
+
+#include <inttypes.h>
+#include <cassert>
+
+namespace chilbert {
+
+// This must be an unsigned integer that is either
+// 32 or 64 bits.  Otherwise, there are places in the
+// code that simply will not work.
+// For speed, this should be the native word size.
+typedef uint64_t FBV_UINT;
+#define FBV_BITS		64
+
+#define FBV0		((FBV_UINT)0)
+#define FBV1		((FBV_UINT)1)
+#define FBV1S		(~FBV0)
+#define FBVN1S(n)	(n==FBV_BITS?FBV1S:(FBV1<<n)-1)
+
+class CFixBitVec
+{
+public:
+
+	// Default constructor.  The bits parameter
+	// is completely ignored, but accepted in order
+	// to look and feel the same as a BigBitVec.
+	CFixBitVec(
+		int iBits = FBV_BITS
+	           )
+	{
+		return;
+	}
+
+	// Copy constructor.
+	CFixBitVec(
+		const CFixBitVec &cFBV
+	           )
+	{
+		m_uiRack = cFBV.m_uiRack;
+	}
+
+	// Returns the current size in bits.
+	int
+	size()
+	{
+		return FBV_BITS;
+	}
+
+	// Zeros the bit-vector.
+	CFixBitVec &
+	reset()
+	{
+		m_uiRack = 0;
+		return (*this);
+	}
+
+	// Truncates the bit-vector to a given precision in
+	// bits (zeroes MSBs without shrinking the vector)
+	CFixBitVec &
+	truncate(
+		int iBits
+	         )
+	{
+		assert( 0 <= iBits && iBits <= FBV_BITS );
+		m_uiRack &= FBVN1S(iBits);
+		return (*this);
+	}
+
+	// Assignment operator.
+	CFixBitVec &
+	operator=(
+		const CFixBitVec &cFBV
+	          )
+	{
+		m_uiRack = cFBV.m_uiRack;
+		return (*this);
+	}
+
+	// Assignment operator.
+	CFixBitVec &
+	operator=(
+		FBV_UINT i
+	          )
+	{
+		m_uiRack = i;
+		return (*this);
+	}
+
+	// Returns the value of the nth bit.
+	bool
+	test(
+		int iIndex
+	       ) const
+	{
+		assert( 0 <= iIndex && iIndex < FBV_BITS );
+		return ( (m_uiRack & (FBV1<<iIndex)) > 0 );
+	}
+
+	// Sets the value of the nth bit.
+	CFixBitVec &
+	set(
+		int iIndex,
+		bool bBit = true
+	       )
+	{
+		assert( 0 <= iIndex && iIndex < FBV_BITS );
+		FBV_UINT m = (FBV1<<iIndex);
+		m_uiRack |= m;
+		if ( !bBit ) m_uiRack ^= m;
+		return (*this);
+	}
+
+	// Toggles the value of the nth bit.
+	CFixBitVec &
+	toggle(
+		int iIndex
+	          )
+	{
+		assert( 0 <= iIndex && iIndex < FBV_BITS );
+		m_uiRack ^= (FBV1<<iIndex);
+		return (*this);
+	}
+
+	// AND operation in place.
+	CFixBitVec &
+	operator&=(
+		const CFixBitVec &cFBV
+	           )
+	{
+		m_uiRack &= cFBV.m_uiRack;
+		return (*this);
+	}
+	CFixBitVec &
+	operator&=(
+		FBV_UINT i
+	           )
+	{
+		m_uiRack &= i;
+		return (*this);
+	}
+
+	// AND operation.
+	CFixBitVec
+	operator&(
+		const CFixBitVec &cFBV
+	          ) const
+	{
+		CFixBitVec t(*this);
+		t &= cFBV;
+		return t;
+	}
+	CFixBitVec
+	operator&(
+		FBV_UINT i
+	          )
+	{
+		CFixBitVec t(*this);
+		t &= i;
+		return t;
+	}
+
+	// OR operation in place.
+	CFixBitVec &
+	operator|=(
+		const CFixBitVec &cFBV
+	           )
+	{
+		m_uiRack |= cFBV.m_uiRack;
+		return (*this);
+	}
+	CFixBitVec &
+	operator|=(
+		FBV_UINT i
+	           )
+	{
+		m_uiRack |= i;
+		return (*this);
+	}
+
+	// OR operation.
+	CFixBitVec
+	operator|(
+		const CFixBitVec &cFBV
+	          ) const
+	{
+		CFixBitVec t(*this);
+		t |= cFBV;
+		return t;
+	}
+	CFixBitVec
+	operator|(
+		FBV_UINT i
+	          )
+	{
+		CFixBitVec t(*this);
+		t |= i;
+		return t;
+	}
+
+
+	// XOR operation in place.
+	CFixBitVec &
+	operator^=(
+		const CFixBitVec &cFBV
+	           )
+	{
+		m_uiRack ^= cFBV.m_uiRack;
+		return (*this);
+	}
+	CFixBitVec &
+	operator^=(
+		FBV_UINT i
+	           )
+	{
+		m_uiRack ^= i;
+		return (*this);
+	}
+
+	// XOR operation.
+	CFixBitVec
+	operator^(
+		const CFixBitVec &cFBV
+	          ) const
+	{
+		CFixBitVec t(*this);
+		t ^= cFBV;
+		return t;
+	}
+	CFixBitVec
+	operator^(
+		FBV_UINT i
+	          )
+	{
+		CFixBitVec t(*this);
+		t ^= i;
+		return t;
+	}
+
+	// Shift left operation, in place.
+	CFixBitVec &
+	operator<<=(
+		int iBits
+	            )
+	{
+		m_uiRack <<= iBits;
+		return (*this);
+	}
+
+	// Shift left operation.
+	CFixBitVec
+	operator<<(
+		int iBits
+	           ) const
+	{
+		CFixBitVec t(*this);
+		t <<= iBits;
+		return t;
+	}
+
+	// Shift right operation, in place.
+	CFixBitVec &
+	operator>>=(
+		int iBits
+	            )
+	{
+		m_uiRack >>= iBits;
+		return (*this);
+	}
+
+	// Shift right operation.
+	CFixBitVec
+	operator>>(
+		int iBits
+	           ) const
+	{
+		CFixBitVec t(*this);
+		t >>= iBits;
+		return t;
+	}
+
+	// Right rotation, in place.
+	CFixBitVec &
+	rotr(
+		int iBits,
+		int iWidth
+	     )
+	{
+		assert( iBits >= 0 );
+		assert( iWidth >  0 );
+		assert( iBits < iWidth );
+		m_uiRack &= FBVN1S(iWidth);
+		m_uiRack = (m_uiRack>>iBits) | (m_uiRack<<(iWidth-iBits));
+		m_uiRack &= FBVN1S(iWidth);
+		return (*this);
+	}
+
+	// Left rotation, in place.
+	CFixBitVec &
+	rotl(
+		int iBits,
+		int iWidth
+	     )
+	{
+		assert( iBits >= 0 );
+		assert( iWidth > 0 );
+		assert( iBits < iWidth );
+		m_uiRack &= FBVN1S(iWidth);
+		m_uiRack = (m_uiRack<<iBits) | (m_uiRack>>(iWidth-iBits));
+		m_uiRack &= FBVN1S(iWidth);
+		return (*this);
+	}
+
+	// Is the bit rack zero valued?
+	bool
+	none() const
+	{
+		return m_uiRack == 0;
+	}
+
+	// Returns the index of the first set bit, numbered from
+	// 1 to n.  0 means there were no set bits.
+	int
+	fsb() const
+	{
+		return chilbert::ffs(m_uiRack);
+	}
+
+
+	// Calculates the Gray Code.
+	CFixBitVec &
+	grayCode()
+	{
+		m_uiRack ^= (m_uiRack>>1);
+		return (*this);
+	}
+
+	// Calculates the Gray Code Inverse
+	CFixBitVec &
+	grayCodeInv()
+	{
+		m_uiRack ^= (m_uiRack>>1);
+		m_uiRack ^= (m_uiRack>>2);
+		m_uiRack ^= (m_uiRack>>4);
+		m_uiRack ^= (m_uiRack>> 8);
+		m_uiRack ^= (m_uiRack>>16);
+#if FBV_BITS == 64
+		m_uiRack ^= (m_uiRack>>32);
+#endif
+		return (*this);
+	}
+
+	// Ones-complements the rack
+	CFixBitVec &
+	flip()
+	{
+		m_uiRack = ~m_uiRack;
+		return (*this);
+	}
+
+	// Returns the first rack.
+	FBV_UINT &
+	rack()
+	{
+		return m_uiRack;
+	}
+	FBV_UINT
+	rack() const
+	{
+		return m_uiRack;
+	}
+
+	// Return a pointer to the racks
+	FBV_UINT *
+	racks()
+	{
+		return &m_uiRack;
+	}
+	const FBV_UINT *
+	racks() const
+	{
+		return &m_uiRack;
+	}
+
+	// Returns the number of racks.
+	int
+	rackCount()
+	{
+		return 1;
+	}
+
+private:
+	static_assert( 8*sizeof(FBV_UINT) == FBV_BITS, "" );
+	static_assert( (sizeof(FBV_UINT) == 4) || (sizeof(FBV_UINT) == 8), "" );
+
+	FBV_UINT	m_uiRack;
+};
+
+} // namespace chilbert
+
+#endif
diff --git a/chilbert/GetBits.hpp b/chilbert/GetBits.hpp
new file mode 100644
index 0000000..48428fc
--- /dev/null
+++ b/chilbert/GetBits.hpp
@@ -0,0 +1,39 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _GETBITS_HPP_
+#define _GETBITS_HPP_
+
+#include "chilbert/Operations.hpp"
+
+namespace chilbert {
+
+template <class H, class I>
+inline void getBits(const H&  h, // bits to read
+                    const int n, // number of bits
+                    const int i, // bit index
+                    I&        w)        // destination
+{
+	for (int j = 0; j < n; j++) {
+		setBit(w, j, testBit(h, i + j));
+	}
+}
+
+} // namespace chilbert
+
+#endif
diff --git a/chilbert/GetLocation.hpp b/chilbert/GetLocation.hpp
new file mode 100644
index 0000000..f41171d
--- /dev/null
+++ b/chilbert/GetLocation.hpp
@@ -0,0 +1,37 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _GETLOCATION_HPP_
+#define _GETLOCATION_HPP_
+
+#include "chilbert/Operations.hpp"
+
+namespace chilbert {
+
+template <class P, class I>
+inline void
+getLocation(const P* const p, const int n, const int i, I& l)
+{
+	for (int j = 0; j < n; ++j) {
+		setBit(l, j, testBit(p[j], i));
+	}
+}
+
+} // namespace chilbert
+
+#endif
diff --git a/chilbert/GrayCodeRank.hpp b/chilbert/GrayCodeRank.hpp
new file mode 100644
index 0000000..13c169c
--- /dev/null
+++ b/chilbert/GrayCodeRank.hpp
@@ -0,0 +1,205 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _GRAYCODERANK_HPP_
+#define _GRAYCODERANK_HPP_
+
+
+#include "chilbert/BigBitVec.hpp"
+
+
+namespace chilbert
+{
+	// This is the bulk of the cost in calculating
+	// a Compact Hilbert Index.  It compresses a previously
+	// calculated index when provided the rotation
+	// at each level of precision.
+	template<class H,class HC>
+	inline
+	void
+	compactIndex(
+		const int *ms,
+		const int *ds,
+		int n,
+		int m,
+		H &h,
+		HC &hc
+		)
+	{
+		int i, j, hr, hcr;
+		FBV_UINT hm, hcm;
+
+		hc.reset();
+
+		hr = hcr = 0;
+		hm = hcm = 1;
+
+		// Run through the levels of precision
+		for ( i = 0; i < m; i++ )
+		{
+			// Run through the dimensions
+			j = ds[i];
+			do
+			{
+				// This dimension contributes a bit?
+				if ( ms[j] > i )
+				{
+					if ( h.racks()[hr] & hm )
+						hc.racks()[hcr] |= hcm;
+					hcm<<=1; if (hcm==0) {hcm=1;++hcr;}
+				}
+
+				if ( ++j == n ) j = 0;
+				hm<<=1; if (hm==0) {hm=1;++hr;}
+			}
+			while ( j != ds[i] );
+		}
+
+		return;
+	}
+
+	template<class I>
+	inline
+	void
+	grayCodeRank(
+		const I &mask,
+		const I &gi,
+		int n,
+		I &r
+		)
+	{
+		r.reset();
+		int i, ir, jr;
+		FBV_UINT im, jm;
+
+		jr = 0; jm = 1;
+		ir = 0; im = 1;
+		for ( i = 0; i < n; ++i )
+		{
+			if ( mask.racks()[ir] & im )
+			{
+				if ( gi.racks()[ir] & im )
+					r.racks()[jr] |= jm;
+				jm<<=1; if ( jm == 0 ) { jm = 1; ++jr; }
+			}
+			
+			im<<=1; if ( im == 0 ) { im = 1; ++ir; }
+		}
+
+		return;
+	}
+
+
+	template<class I>
+	inline
+	void
+	grayCodeRankInv(
+		const I &mask,
+		const I &ptrn,
+		const I &r,
+		int n,
+		int b,
+		I &g,
+		I &gi
+		)
+	{
+		g.reset();
+		gi.reset();
+
+		int i, j, ir, jr;
+		FBV_UINT im, jm;
+
+		i = n-1;
+		BBV_MODSPLIT(ir,im,i);
+		im = (FBV1<<im);
+
+		j = b-1;
+		BBV_MODSPLIT(jr,jm,j);
+		jm = (FBV1<<jm);
+
+		FBV_UINT gi0, gi1, g0;
+		gi1 = gi0 = g0 = 0;
+
+		for ( ; i >= 0; --i )
+		{
+			// Unconstrained bit?
+			if ( mask.racks()[ir] & im )
+			{
+				gi1 = gi0;
+				gi0 = (r.racks()[jr] & jm) > 0;
+				g0 = gi0 ^ gi1;
+				if ( gi0 ) gi.racks()[ir] |= im;
+				if ( g0 ) g.racks()[ir] |= im;
+				jm>>=1; if ( jm == 0 ) { jm=((FBV_UINT)1)<<(FBV_BITS-1); --jr; }
+			}
+			else
+			{
+				g0 = (ptrn.racks()[ir] & im) > 0;
+				gi1 = gi0;
+				gi0 = g0 ^ gi1;
+				if ( gi0 ) gi.racks()[ir] |= im;
+				if ( g0 ) g.racks()[ir] |= im;
+			}
+
+			im>>=1; if ( im == 0 ) { im=((FBV_UINT)1)<<(FBV_BITS-1); --ir; }
+		}
+
+		return;
+	}
+
+	
+	template <class I>
+	inline
+	void
+	extractMask(
+		const int *ms,
+		int n,
+		int d,
+		int i,
+		I &mask,
+		int &b
+		)
+	{
+		assert( 0 <= d && d < n );
+		int j, jr;
+		FBV_UINT jm;
+
+		mask.reset();
+		b = 0;
+		
+		jm = 1; jr = 0;
+		j = d; // #D j = (d==n-1) ? 0 : d+1;
+		do
+		{
+			if ( ms[j] > i )
+			{
+				mask.racks()[jr] |= jm;
+				++b;
+			}
+		
+			jm <<= 1; if ( jm == 0 ) { jm=1; ++jr; }
+			if ( ++j == n ) j = 0;
+		}
+		while ( j != d );
+
+		return;
+	}
+};
+
+
+#endif
diff --git a/chilbert/Hilbert.hpp b/chilbert/Hilbert.hpp
new file mode 100644
index 0000000..b6d41b0
--- /dev/null
+++ b/chilbert/Hilbert.hpp
@@ -0,0 +1,141 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _HILBERT_HPP_
+#define _HILBERT_HPP_
+
+
+#include "chilbert/FixBitVec.hpp"
+#include "chilbert/BigBitVec.hpp"
+#include "chilbert/Algorithm.hpp"
+
+// Description of parameters:
+//
+// FOR REGULAR HILBERT INDICES
+//
+// CFixBitVec/CBigBitVec *p
+// Pointer to array of non-negative coordinate values.
+//
+// int m
+// Precision of all coordinate values (number of bits required to
+// represent the largest possible coordinate value).
+//
+// int n
+// Number of dimensions (size of the array *p).
+//
+// CFixBitVec/CBigBitVec &h
+// Hilbert index of maximum precision m*n.
+//
+// int *ms
+// Array of precision values, one per dimension.
+//
+// FOR COMPACT HILBERT INDICES
+//
+// CFixBitVec/CBigBitVec &hc
+// Compact Hilbert index of maximum precision M.
+//
+// int M
+// Net precision value, corresponding to the size of the compact
+// Hilbert code.  If not provided, defaults to zero and will be calculated
+// by the function (sum_i { ms[i] }).
+//
+// int m
+// Largest precision value (max_i { ms[i] }).  If not provided, defaults
+// to zero and will be calculated by the function,
+
+
+namespace chilbert
+{
+
+// fix -> fix
+
+inline void coordsToIndex( const CFixBitVec *p, int m, int n, CFixBitVec &h )
+{
+	coordsToIndex<CFixBitVec,CFixBitVec>(p,m,n,h);
+}
+
+inline void indexToCoords( CFixBitVec *p, int m, int n, const CFixBitVec &h )
+{
+	indexToCoords<CFixBitVec,CFixBitVec>(p,m,n,h);
+}
+
+inline void coordsToCompactIndex( const CFixBitVec *p, const int *ms, int n,
+                                  CFixBitVec &hc, int M, int m )
+{
+	coordsToCompactIndex<CFixBitVec,CFixBitVec>(p,ms,n,hc,M,m);
+}
+
+inline void compactIndexToCoords( CFixBitVec *p, const int *ms, int n,
+                                  const CFixBitVec &hc, int M, int m )
+{
+	compactIndexToCoords<CFixBitVec,CFixBitVec>(p,ms,n,hc,M,m);
+}
+
+// fix -> big
+
+inline void coordsToIndex( const CFixBitVec *p, int m, int n, CBigBitVec &h )
+{
+	coordsToIndex<CFixBitVec,CBigBitVec>(p,m,n,h);
+}
+
+inline void indexToCoords( CFixBitVec *p, int m, int n, const CBigBitVec &h )
+{
+	indexToCoords<CFixBitVec,CBigBitVec>(p,m,n,h);
+}
+
+inline void coordsToCompactIndex( const CFixBitVec *p, const int *ms, int n,
+                                  CBigBitVec &hc, int M, int m )
+{
+	coordsToCompactIndex<CFixBitVec,CBigBitVec>(p,ms,n,hc,M,m);
+}
+
+inline void compactIndexToCoords( CFixBitVec *p, const int *ms, int n,
+                                  const CBigBitVec &hc, int M, int m )
+{
+	compactIndexToCoords<CFixBitVec,CBigBitVec>(p,ms,n,hc,M,m);
+}
+
+// big -> big
+
+inline void coordsToIndex( const CBigBitVec *p, int m, int n, CBigBitVec &h )
+{
+	coordsToIndex<CBigBitVec,CBigBitVec>(p,m,n,h);
+}
+
+inline void indexToCoords( CBigBitVec *p, int m, int n, const CBigBitVec &h )
+{
+	indexToCoords<CBigBitVec,CBigBitVec>(p,m,n,h);
+}
+
+inline void coordsToCompactIndex( const CBigBitVec *p, const int *ms, int n,
+                                  CBigBitVec &hc, int M, int m )
+{
+	coordsToCompactIndex<CBigBitVec,CBigBitVec>(p,ms,n,hc,M,m);
+}
+
+inline void compactIndexToCoords( CBigBitVec *p, const int *ms, int n,
+                                  const CBigBitVec &hc, int M, int m )
+{
+	compactIndexToCoords<CBigBitVec,CBigBitVec>(p,ms,n,hc,M,m);
+}
+
+};
+
+
+#endif
+
diff --git a/chilbert/Operations.hpp b/chilbert/Operations.hpp
new file mode 100644
index 0000000..01f9197
--- /dev/null
+++ b/chilbert/Operations.hpp
@@ -0,0 +1,92 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _OPERATIONS_HPP_
+#define _OPERATIONS_HPP_
+
+#include <cassert>
+#include <climits>
+#include <cstddef>
+#include <type_traits>
+
+namespace chilbert {
+
+/// IntegralIndex<T> only exists if T is integral
+template <typename T>
+using IntegralIndex = std::enable_if_t<std::is_integral<T>::value, int>;
+
+/// BitsetIndex<T> only exists if T is not integral (must be a bitset)
+template <typename T>
+using BitsetIndex = std::enable_if_t<!std::is_integral<T>::value, int>;
+
+/// Return the `index`th bit in `field`
+template <typename T>
+bool
+testBit(const T& field, const IntegralIndex<T> index)
+{
+	assert(size_t(index) < sizeof(field) * CHAR_BIT);
+	return field & (((T)1) << index);
+}
+
+/// Return the `index`th bit in `field`
+template <typename T>
+bool
+testBit(const T& field, const BitsetIndex<T> index)
+{
+	return field.test(index);
+}
+
+/// Set the `index`th bit in `field` to `value`
+template <typename T>
+void
+setBit(T& field, const IntegralIndex<T> index, const bool value)
+{
+	assert(size_t(index) < sizeof(field) * CHAR_BIT);
+	field ^= (-value ^ field) & ((T)1 << index);
+	assert(testBit(field, index) == value);
+}
+
+/// Set the `index`th bit in `field` to `value`
+template <typename T>
+void
+setBit(T& field, const BitsetIndex<T> index, const bool value)
+{
+	field.set(index, value);
+}
+
+/// Return 1 + the index of the least significant 1-bit of `field`, or zero
+template <typename T>
+int ffs(const T field);
+
+template <>
+int
+ffs<unsigned long>(const unsigned long field)
+{
+	return __builtin_ffsl(field);
+}
+
+template <>
+int
+ffs<unsigned long long>(const unsigned long long field)
+{
+	return __builtin_ffsll(field);
+}
+
+} // namespace chilbert
+
+#endif
diff --git a/chilbert/SetBits.hpp b/chilbert/SetBits.hpp
new file mode 100644
index 0000000..adca0a0
--- /dev/null
+++ b/chilbert/SetBits.hpp
@@ -0,0 +1,39 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _SETBITS_HPP_
+#define _SETBITS_HPP_
+
+#include "chilbert/Operations.hpp"
+
+namespace chilbert {
+
+template <class H, class I>
+inline void setBits(H&        h, // destination
+                    const int n, // number of bits
+                    const int i, // bit position
+                    const I&  w)  // bits to place
+{
+	for (int j = 0; j < n; j++) {
+		setBit(h, i + j, testBit(w, j));
+	}
+}
+
+} // namespace chilbert
+
+#endif
diff --git a/chilbert/SetLocation.hpp b/chilbert/SetLocation.hpp
new file mode 100644
index 0000000..d8dec11
--- /dev/null
+++ b/chilbert/SetLocation.hpp
@@ -0,0 +1,37 @@
+/*
+  Copyright (C) 2018 David Robillard <d@drobilla.net>
+  Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+
+  This program is free software: you can redistribute it and/or modify it under
+  the terms of the GNU General Public License as published by the Free Software
+  Foundation, either version 2 of the License, or (at your option) any later
+  version.
+
+  This program 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 General Public License for more
+  details.
+
+  You should have received a copy of the GNU General Public License along with
+  this program.  If not, see <https://www.gnu.org/licenses/>.
+*/
+
+#ifndef _SETLOCATION_HPP_
+#define _SETLOCATION_HPP_
+
+#include "chilbert/Operations.hpp"
+
+namespace chilbert {
+
+template <class P, class I>
+inline void
+setLocation(P* const p, const int n, const int i, const I& l)
+{
+	for (int j = 0; j < n; j++) {
+		setBit(p[j], i, testBit(l, j));
+	}
+}
+
+} // namespace chilbert
+
+#endif