Initial import

10 files changed, 2758 insertions, 0 deletions

diff --git a/Hilbert/Algorithm.hpp b/Hilbert/Algorithm.hpp
new file mode 100644
index 0000000..343774c
--- /dev/null
+++ b/Hilbert/Algorithm.hpp
@@ -0,0 +1,482 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _ALGORITHM_HPP_
+#define _ALGORITHM_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/BigBitVec.hpp>
+#include <Hilbert/GetLocation.hpp>
+#include <Hilbert/SetLocation.hpp>
+#include <Hilbert/GetBits.hpp>
+#include <Hilbert/SetBits.hpp>
+#include <Hilbert/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 Hilbert
+{
+	// 'Transforms' a point.
+	template<class I>
+	H_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>
+	H_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>
+	H_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.toggleBit( 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.toggleBit( d == 0 ? n-1 : d-1 ); //#D d );
+
+		return;
+	}
+
+	// Update for method 2 (GrayCodeInv out of loop)
+	template<class I>
+	H_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.toggleBit( 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>
+	H_INLINE
+	void
+	_coordsToIndex(
+		const P *p,
+		int m,
+		int n,
+		H &h,
+		int *ds = NULL // #HACK
+		)
+	{
+		I e(n), l(n), t(n), w(n);
+		int d, i;
+		int ho = m*n;
+
+		// Initialize
+		e.zero();
+		d = D0;
+		l.zero();
+		h.zero();
+
+		int r, b;
+		BBV_MODSPLIT(r,b,n-1);
+		FBV_UINT bm = (1<<b);
+
+		// 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.racks()[r] ^= bm;
+
+			// 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>
+	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>
+	H_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.zero();
+		d = D0;
+		l.zero();
+		for ( j = 0; j < n; j++ )
+			p[j].zero();
+
+		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>
+	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>
+	H_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>
+	H_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>
+	H_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.zero();
+		d = D0;
+		l.zero();
+		for ( j = 0; j < n; j++ )
+			p[j].zero();
+		
+		// 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.zero(); // 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>
+	H_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/Hilbert/BigBitVec.hpp b/Hilbert/BigBitVec.hpp
new file mode 100644
index 0000000..2bb6db7
--- /dev/null
+++ b/Hilbert/BigBitVec.hpp
@@ -0,0 +1,928 @@
+/*
+ * Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+ * Copyright (C) 2016 David Robillard <d@drobilla.net>
+ *
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _BIGBITVEC_HPP_
+#define _BIGBITVEC_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/FixBitVec.hpp>
+
+#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; }
+
+class CBigBitVec
+{
+public:
+
+	static EBitVecType
+	type()
+	{
+		return eBig;
+	}
+
+
+	// 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 = NULL;
+		}
+			
+		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 = NULL;
+	}
+
+	// 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
+	getSize() const
+	{
+		return m_iRacks*FBV_BITS;
+	}
+
+
+	// Resize function.	Returns the number of bits
+	// we can accomodate after resizing.
+	CBigBitVec &
+	setSize(
+		int iBits
+	        )
+	{
+		int i;
+	
+		// How many racks do we need?
+		int iRacks = FBVS_NEEDED(iBits);
+
+		// Same size?
+		if ( iRacks == m_iRacks )
+			return (*this);
+	
+		// Allocate new racks.
+		CFixBitVec *pcRacks = new CFixBitVec[iRacks];
+
+		// Copy over the old values.
+		/*for ( i = 0; i < BBV_MIN(iRacks,m_iRacks); i++ )
+		  pcRacks[i] = m_pcRacks[i];*/
+		i = BBV_MIN(iRacks,m_iRacks);
+		if (m_pcRacks) {
+			memcpy( static_cast<void*>(pcRacks),
+			        static_cast<const void*>(m_pcRacks),
+			        sizeof(CFixBitVec)*i );
+		}
+
+		// zero the new values
+		/*for ( ; i < iRacks; i++ )
+		  pcRacks[i].zero();*/
+		if ( iRacks > i )
+			memset( static_cast<void*>(pcRacks + i), 0,
+			        sizeof(CFixBitVec)*(iRacks-i) );
+
+		// Release the old stuff.
+		delete [] m_pcRacks;
+
+		// Replace old with new.
+		m_iRacks = iRacks;
+		m_pcRacks = pcRacks;
+	
+		return (*this);
+	}
+
+
+	// zeros the bit-vector.
+	CBigBitVec &
+	zero()
+	{
+	
+		/*int i;
+		  for ( i = 0; i < m_iRacks; i++ )
+		  m_pcRacks[i].zero();*/
+		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 <= getSize() );
+		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].zero();
+
+		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].zero();*/
+			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 = NULL;
+		}
+
+		return (*this);
+	}
+	CBigBitVec &
+	operator=(
+		const CFixBitVec &cFBV
+	          )
+	{
+		m_pcRacks[0] = cFBV;
+		/*int i;
+		  for ( i = 1; i < m_iRacks; i++ )
+		  m_pcRacks[i].zero();*/
+		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].zero();*/
+		memset( static_cast<void*>(m_pcRacks+1),
+		        0, sizeof(CFixBitVec)*(m_iRacks-1) );
+		return (*this);
+	}
+
+	// Returns the value of the nth bit.
+	bool
+	getBit(
+		int iIndex
+	       ) const
+	{
+		assert( iIndex >= 0 && iIndex < getSize() );
+		int r, b;
+		BBV_MODSPLIT(r,b,iIndex);
+		return m_pcRacks[r].getBit(b);
+	}
+
+
+	// Sets the value of the nth bit.
+	CBigBitVec &
+	setBit(
+		int iIndex,
+		bool bBit
+	       )
+	{
+		assert( iIndex >= 0 && iIndex < getSize() );
+		int r, b;
+		BBV_MODSPLIT(r,b,iIndex);
+		m_pcRacks[r].setBit(b,bBit);
+		return (*this);
+	}
+
+
+	// Toggles the value of the nth bit.
+	CBigBitVec &
+	toggleBit(
+		int iIndex
+	          )
+	{
+		assert( iIndex >= 0 && iIndex < getSize() );
+		int r, b;
+		BBV_MODSPLIT(r,b,iIndex);
+		m_pcRacks[r].toggleBit(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 )
+		{
+			zero();
+			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].zero();
+		}
+
+		// 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 )
+		{
+			zero();
+			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].zero();
+		}
+
+		// 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 = getSize();
+
+		// Modulo the number of bits.
+		//FBVMOD(iBits,iWidth);
+		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);
+	}
+
+
+	// Right rotation.
+	CBigBitVec
+	rotrCopy(
+		int iBits,
+		int iWidth
+	         ) const
+	{
+		CBigBitVec t( *this );
+		t.rotr(iBits,iWidth);
+		return t;
+	}
+
+
+	// Left rotation, in place.
+	CBigBitVec &
+	rotl(
+		int iBits,
+		int iWidth
+	     )
+	{
+		assert( iBits >= 0 );
+
+		// Fill in the width, if necessary.
+		if ( iWidth <= 0 )
+			iWidth = getSize();
+
+		// Modulo the number of bits.
+		//FBVMOD(iBits,iWidth);
+		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.
+	CBigBitVec
+	rotlCopy(
+		int iBits,
+		int iWidth
+	         ) const
+	{
+		CBigBitVec t( *this );
+		t.rotl(iBits,iWidth);
+		return t;
+	}
+
+
+	// Returns true if the rack is zero valued.
+	bool
+	isZero() const
+	{
+		int i;
+		for ( i = 0; i < m_iRacks; i++ )
+			if ( !m_pcRacks[i].isZero() ) return false;
+		return true;
+	}
+
+
+	// Returns the number of trailing set bits.
+	int
+	tsb() const
+	{
+		int c, i, j;
+		c = 0;
+		for ( i = 0; i < m_iRacks; i++ )
+		{
+			j = m_pcRacks[i].tsb();
+			c += j;
+			if ( j < FBV_BITS )
+				break;
+		}
+		return c;
+	}
+
+	// OB:
+	// Returns the index of the most significant bit (numbered
+	// 1 to n)
+	int
+	msb() const
+	{
+		int c, i, j = 0;
+		c = FBV_BITS * m_iRacks;
+		for ( i = m_iRacks - 1; i >= 0 && !j; i-- )
+		{
+			j = m_pcRacks[i].msb();
+			if (j)
+				return c - (FBV_BITS - j);
+			else
+				c -= FBV_BITS;
+		}
+		return 0;
+	}
+
+	// Returns the index of the first set bit.
+	// (numbered 1 to n, with 0 meaning no bits were set)
+	int
+	fsb() const
+	{
+		int c, i, j;
+		c = 0;
+		for ( i = 0; i < m_iRacks; i++ )
+		{
+			j = m_pcRacks[i].fsb();
+			if ( j < FBV_BITS )
+				return c + j;
+			else
+				c += FBV_BITS;
+		}
+		return 0;
+	}
+
+
+	// Prefix decrement.	Returns true if there
+	// was a carry, false otherwise.
+	bool
+	operator--()
+	{
+		int i = 0;
+		bool b;
+		while ( i < m_iRacks && (b = --m_pcRacks[i]) ) i++;
+	
+		return b;
+	}
+
+
+	// 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].complement();
+			s = m_pcRacks[i].rack() & 1;
+		}
+		return (*this);
+	}
+
+
+	// Complement
+	CBigBitVec &
+	complement()
+	{
+		int i;
+		for ( i = 0; i < m_iRacks; i++ )
+			m_pcRacks[i].complement();
+		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;
+};
+
+#endif
diff --git a/Hilbert/Common.hpp b/Hilbert/Common.hpp
new file mode 100644
index 0000000..b297162
--- /dev/null
+++ b/Hilbert/Common.hpp
@@ -0,0 +1,26 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _COMMON_HPP_
+#define _COMMON_HPP_
+
+#include <assert.h>
+
+#define H_INLINE	inline
+
+#endif
diff --git a/Hilbert/FixBitVec.hpp b/Hilbert/FixBitVec.hpp
new file mode 100644
index 0000000..84015cf
--- /dev/null
+++ b/Hilbert/FixBitVec.hpp
@@ -0,0 +1,538 @@
+/*
+ * Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+ * Copyright (C) 2015 David Robillard <d@drobilla.net>
+ *
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _FIXBITVEC_HPP_
+#define _FIXBITVEC_HPP_
+
+
+#include <inttypes.h>
+#include <Hilbert/Common.hpp>
+
+
+// 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)
+#define FBVMOD(i,m)	if((i)>=(m))(i)-=(m)*((i)/(m));
+
+#define COMPILE_TIME_ASSERT(pred) switch(0){case 0:case pred:;}
+
+enum EBitVecType
+{
+	eFix,
+	eBig
+};
+
+
+class CFixBitVec
+{
+public:
+
+	static EBitVecType
+	type()
+	{
+		return eFix;
+	}
+
+	// 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
+	getSize()
+	{
+		return FBV_BITS;
+	}
+
+	// Sets the size.  This is a dummy
+	// function just for BigBitVec compatibility.
+	CFixBitVec &
+	setSize(
+		int iBits
+	        )
+	{
+		return (*this);
+	}
+
+	// Zeros the bit-vector.
+	CFixBitVec &
+	zero()
+	{
+		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
+	getBit(
+		int iIndex
+	       ) const
+	{
+		assert( 0 <= iIndex && iIndex < FBV_BITS );
+		return ( (m_uiRack & (FBV1<<iIndex)) > 0 );
+	}
+
+	// Sets the value of the nth bit.
+	CFixBitVec &
+	setBit(
+		int iIndex,
+		bool bBit
+	       )
+	{
+		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 &
+	toggleBit(
+		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 );//#D FBVMOD(iBits,iWidth);
+		m_uiRack &= FBVN1S(iWidth);
+		m_uiRack = (m_uiRack>>iBits) | (m_uiRack<<(iWidth-iBits));
+		m_uiRack &= FBVN1S(iWidth);
+		return (*this);
+	}
+
+	// Right rotation.
+	CFixBitVec
+	rotrCopy(
+		int iBits,
+		int iWidth
+	         ) const
+	{
+		CFixBitVec t(*this);
+		t.rotr(iBits,iWidth);
+		return t;
+	}
+
+	// Left rotation, in place.
+	CFixBitVec &
+	rotl(
+		int iBits,
+		int iWidth
+	     )
+	{
+		assert( iBits >= 0 );
+		assert( iWidth > 0 );
+		assert( iBits < iWidth );//#D FBVMOD(iBits,iWidth);
+		m_uiRack &= FBVN1S(iWidth);
+		m_uiRack = (m_uiRack<<iBits) | (m_uiRack>>(iWidth-iBits));
+		m_uiRack &= FBVN1S(iWidth);
+		return (*this);
+	}
+
+	// Left rotation.
+	CFixBitVec
+	rotlCopy(
+		int iBits,
+		int iWidth
+	         ) const
+	{
+		CFixBitVec t(*this);
+		t.rotl(iBits,iWidth);
+		return t;
+	}
+
+	// Is the bit rack zero valued?
+	bool
+	isZero() const
+	{
+		return m_uiRack == 0;
+	}
+
+	// Returns the number of trailing set bits.
+	int
+	tsb() const
+	{
+		FBV_UINT i = m_uiRack;
+		int c = 0;
+
+#if FBV_BITS == 64
+		if ( i == FBV1S ) return 64;
+		if ( (i&FBVN1S(32)) == FBVN1S(32) ) { i>>=32; c^=32; }
+#elif FBV_BITS == 32
+		if ( i == FBV1S ) return 32;
+#endif
+		if ( (i&FBVN1S(16)) == FBVN1S(16) ) { i>>=16; c^=16; }
+		if ( (i&FBVN1S( 8)) == FBVN1S( 8) ) { i>>= 8; c^= 8; }
+		if ( (i&FBVN1S( 4)) == FBVN1S( 4) ) { i>>= 4; c^= 4; }
+		if ( (i&FBVN1S( 2)) == FBVN1S( 2) ) { i>>= 2; c^= 2; }
+		if ( (i&FBVN1S( 1)) == FBVN1S( 1) ) { i>>= 1; c^= 1; }
+
+		return c;	
+	}
+
+	// Returns the index of the most significant bit
+	int
+	msb() const
+	{
+		FBV_UINT i = m_uiRack;
+		int c = 0;
+#if FBV_BITS == 64
+		if ( i == FBV0 ) return 0;
+		if ( i & (FBVN1S(32) << 32) ) { i >>= 32; c ^= 32; }
+#elif FBV_BITS == 32
+		if ( i == FBV0 ) return 0;
+#endif
+		if ( i & (FBVN1S(16) << 16) ) { i >>= 16; c ^= 16; }
+		if ( i & (FBVN1S( 8) <<  8) ) { i >>=  8; c ^=  8; }
+		if ( i & (FBVN1S( 4) <<  4) ) { i >>=  4; c ^=  4; }
+		if ( i & (FBVN1S( 2) <<  2) ) { i >>=  2; c ^=  2; }
+		if ( i & (FBVN1S( 1) <<  1) ) { i >>=  1; c ^=  1; }
+		return ++c;
+	}
+                
+	// Returns the index of the first set bit, numbered from
+	// 1 to n.  0 means there were no set bits.
+	int
+	fsb() const
+	{
+		FBV_UINT i = m_uiRack;
+		int c = 0;
+
+#if FBV_BITS == 64
+		if ( i == FBV0 ) return 0;
+		if ( (i&FBVN1S(32)) == FBV0 ) { i>>=32; c^=32; }
+#elif FBV_BITS == 32
+		if ( i == FBV0 ) return 0;
+#endif
+		if ( (i&FBVN1S(16)) == FBV0 ) { i>>=16; c^=16; }
+		if ( (i&FBVN1S( 8)) == FBV0 ) { i>>= 8; c^= 8; }
+		if ( (i&FBVN1S( 4)) == FBV0 ) { i>>= 4; c^= 4; }
+		if ( (i&FBVN1S( 2)) == FBV0 ) { i>>= 2; c^= 2; }
+		if ( (i&FBVN1S( 1)) == FBV0 ) { i>>= 1; c^= 1; }
+
+		return ++c;
+	}
+
+	// Prefix decrement.  Returns true if a carry
+	// was generated.
+	bool
+	operator--()
+	{
+		return ( m_uiRack-- == 0 );
+	}
+
+	// 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 &
+	complement()
+	{
+		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 void
+	compileTimeAssertions()
+	{
+		COMPILE_TIME_ASSERT( 8*sizeof(FBV_UINT) == FBV_BITS );
+		COMPILE_TIME_ASSERT( (sizeof(FBV_UINT) == 4) || (sizeof(FBV_UINT) == 8) );
+	}
+	
+	FBV_UINT	m_uiRack;
+};
+
+
+#endif
diff --git a/Hilbert/GetBits.hpp b/Hilbert/GetBits.hpp
new file mode 100644
index 0000000..119ace9
--- /dev/null
+++ b/Hilbert/GetBits.hpp
@@ -0,0 +1,93 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _GETBITS_HPP_
+#define _GETBITS_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/BigBitVec.hpp>
+
+
+namespace Hilbert
+{
+	template <class H,class I>
+	H_INLINE
+	void
+	getBits(
+		const H &h,	// bits to read
+		int n,			// number of bits
+		int i,			// bit index
+		I &w				// destination
+		)
+	{
+		// This is terribly inefficient.
+		int j;
+		for ( j = 0; j < n; j++ )
+			w.setBit(j,h.getBit(i+j));
+		return;
+	}
+
+	
+	// <CBigBitVec,CBigBitVec>
+	// #TODO
+
+	// <CBigBitVec,CFixBitVec>
+	template<>
+	H_INLINE
+	void
+	getBits(
+		const CBigBitVec &h,
+		int n,
+		int i,
+		CFixBitVec &w
+		)
+	{
+		int ir, ib, t;
+		BBV_MODSPLIT(ir,ib,i);
+		w.rack() = h.racks()[ir] >> ib;
+		t = FBV_BITS - ib;
+		if ( t < n )
+		{
+			w.rack() |= h.racks()[ir+1] >> (FBV_BITS-n);
+		}
+		w.truncate(n);
+		return;
+	}
+	
+	
+	// <CFixBitVec,CFixBitVec>
+	template<>
+	H_INLINE
+	void
+	getBits(
+		const CFixBitVec &h,
+		int n,
+		int i,
+		CFixBitVec &w
+		)
+	{
+		w = h >> i;
+		w.truncate(n);
+		return;
+	}
+
+};
+
+
+#endif
diff --git a/Hilbert/GetLocation.hpp b/Hilbert/GetLocation.hpp
new file mode 100644
index 0000000..99ea7ba
--- /dev/null
+++ b/Hilbert/GetLocation.hpp
@@ -0,0 +1,114 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _GETLOCATION_HPP_
+#define _GETLOCATION_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/BigBitVec.hpp>
+
+
+namespace Hilbert
+{
+	
+	template<class P,class I>
+	H_INLINE
+	void
+	_getLocation(
+		const P *p,
+		int jo,
+		int jn,
+		int ir,
+		FBV_UINT im,
+		FBV_UINT &l
+		)
+	{
+		l = 0;
+    switch ( jn )
+    {
+#define GETLOC_CASE(i)		case ((i)+1): if (p[jo+(i)].racks()[ir]&im) l|=(FBV1<<(i))
+#define GETLOC_CASE2(i) \
+			GETLOC_CASE(i+1); \
+			GETLOC_CASE(i)
+#define GETLOC_CASE4(i) \
+			GETLOC_CASE2(i+2); \
+			GETLOC_CASE2(i)
+#define GETLOC_CASE8(i) \
+			GETLOC_CASE4(i+4); \
+			GETLOC_CASE4(i)
+#define GETLOC_CASE16(i) \
+			GETLOC_CASE8(i+8); \
+			GETLOC_CASE8(i)
+#define GETLOC_CASE32(i) \
+			GETLOC_CASE16(i+16); \
+			GETLOC_CASE16(i)
+#if FBV_BITS == 64
+			GETLOC_CASE32(32);
+#endif
+			GETLOC_CASE32(0);
+		}
+		return;
+	}
+	
+	
+	template<class P,class I>
+	H_INLINE
+	void
+	getLocation(
+		const P *p,
+		int n,
+		int i,
+		I	&l
+		)
+	{
+		/*int j;
+		for ( j = n-1; j >= 0; --j )
+			l.setBit(j,p[j].getBit(i));
+		return;*/
+		
+		int j, jo, ir;
+		FBV_UINT im;
+		
+		if ( P::type() == eBig )
+		{
+			ir = i / FBV_BITS;
+			im = FBV1 << (i-ir*FBV_BITS);
+		}
+		else
+		{
+			ir = 0;
+			im = FBV1 << i;
+		}
+		
+		j = 0;
+		jo = 0;
+		if ( I::type() == eBig )
+		{
+			for ( ; j < l.rackCount()-1; ++j, jo += FBV_BITS )
+				_getLocation<P,I>(p,jo,FBV_BITS,ir,im,l.racks()[j]);
+		}
+		_getLocation<P,I>(p,jo,n-jo,ir,im,l.racks()[j]);
+		
+		return;
+	}
+	
+};
+
+
+#endif
diff --git a/Hilbert/GrayCodeRank.hpp b/Hilbert/GrayCodeRank.hpp
new file mode 100644
index 0000000..33ea9f4
--- /dev/null
+++ b/Hilbert/GrayCodeRank.hpp
@@ -0,0 +1,246 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _GRAYCODERANK_HPP_
+#define _GRAYCODERANK_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/BigBitVec.hpp>
+
+
+namespace Hilbert
+{
+	// 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>
+	H_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.zero();
+
+		hr = hcr = 0;
+		hm = hcm = 1;
+
+		// Run through the levels of precision
+		for ( i = 0; i < m; i++ )
+		{
+			// #D
+			// int k = ds[i] + 1;
+			// if ( k == n ) k = 0;
+			// j = k;
+
+/*#define COMP	{ \
+	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;} \
+	}
+#define COMP1(a)	case a: COMP;
+#define COMP2(a)	COMP1(a+1); \
+		COMP1(a);
+#define COMP4(a)	COMP2(a+2); \
+		COMP2(a);
+#define COMP8(a)	COMP4(a+4); \
+		COMP4(a);
+#define COMP16(a)	COMP8(a+8); \
+		COMP8(a);
+#define COMP32(a)	COMP16(a+16); \
+		COMP16(a);
+
+			// This complicated mess only buys a marginal performance increase.
+			int k, kd;
+			k = n;
+			j = ds[i];
+			while ( k )
+			{
+				kd = k;
+				switch ( kd )
+				{
+					default: kd = 32;
+					COMP32(1);
+				}
+				k -= kd;
+			}*/
+			
+			// 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>
+	H_INLINE
+	void
+	grayCodeRank(
+		const I &mask,
+		const I &gi,
+		int n,
+		I &r
+		)
+	{
+		r.zero();
+		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>
+	H_INLINE
+	void
+	grayCodeRankInv(
+		const I &mask,
+		const I &ptrn,
+		const I &r,
+		int n,
+		int b,
+		I &g,
+		I &gi
+		)
+	{
+		g.zero();
+		gi.zero();
+
+		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>
+	H_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.zero();
+		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/Hilbert/Hilbert.hpp b/Hilbert/Hilbert.hpp
new file mode 100644
index 0000000..bd39962
--- /dev/null
+++ b/Hilbert/Hilbert.hpp
@@ -0,0 +1,142 @@
+/*
+ * Copyright (C) 2006-2007 Chris Hamilton <chamilton@cs.dal.ca>
+ * Copyright (C) 2015 David Robillard <d@drobilla.net>
+ *
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _HILBERT_HPP_
+#define _HILBERT_HPP_
+
+
+#include <Hilbert/FixBitVec.hpp>
+#include <Hilbert/BigBitVec.hpp>
+#include <Hilbert/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 Hilbert
+{
+
+// 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/Hilbert/SetBits.hpp b/Hilbert/SetBits.hpp
new file mode 100644
index 0000000..87f19cb
--- /dev/null
+++ b/Hilbert/SetBits.hpp
@@ -0,0 +1,93 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _SETBITS_HPP_
+#define _SETBITS_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/BigBitVec.hpp>
+
+
+namespace Hilbert
+{
+
+	template<class H,class I>
+	H_INLINE
+	void
+	setBits(
+		H &h,				// destination
+		int n,			// number of bits
+		int i,			// bit position
+		const I &w	// bits to place
+		)
+	{
+		// This is terribly inefficient.
+		int j;
+		for ( j = 0; j < n; j++ )
+			h.setBit(i+j,w.getBit(j));
+	}
+
+
+	// <CBigBitVec,CBigBitVec>
+	// #TODO
+
+
+	// <CBigBitVec,CFixBitVec>
+	template<>
+	H_INLINE
+	void
+	setBits(
+		CBigBitVec &h,
+		int n,
+		int i,
+		const CFixBitVec &w
+		)
+	{
+		int ir, ib, t;
+		BBV_MODSPLIT(ir,ib,i);
+		h.racks()[ir] |= w.rack() << ib;
+		t = ib+n;
+		if ( t > FBV_BITS )
+		{
+			t -= FBV_BITS;
+			h.racks()[ir+1] |= w.rack() >> t;
+		}
+		return;
+	}
+	
+	
+	// <CFixBitVec,CFixBitVec>
+	template<>
+	H_INLINE
+	void
+	setBits(
+		CFixBitVec &h,
+		int n,
+		int i,
+		const CFixBitVec &w
+		)
+	{
+		h |= w << i;
+		return;
+	}
+
+};
+
+
+#endif
diff --git a/Hilbert/SetLocation.hpp b/Hilbert/SetLocation.hpp
new file mode 100644
index 0000000..d7b5099
--- /dev/null
+++ b/Hilbert/SetLocation.hpp
@@ -0,0 +1,96 @@
+/*
+ * 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, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
+#ifndef _SETLOCATION_HPP_
+#define _SETLOCATION_HPP_
+
+
+#include <Hilbert/Common.hpp>
+#include <Hilbert/BigBitVec.hpp>
+
+
+namespace Hilbert
+{
+	template <class P,class I>
+	H_INLINE
+	void
+	setLocation(
+		P *p,
+		int n,
+		int i,
+		const I &l
+		)
+	{
+		// Easy to understand implementation
+		/*int j;
+		for ( j = 0; j < n; j++ )
+			p[j].setBit(i,l.getBit(j));
+		return;*/
+
+		// Much faster loop-unrolled implementation.
+		int ir, ib;
+		FBV_UINT im;
+		BBV_MODSPLIT(ir,ib,i);
+		im = (FBV1<<ib);
+
+#define SETBIT p->racks()[ir]|=im; if ((*lr&jm)==0) p->racks()[ir]^=im; jm<<=1; ++p;
+#define SETBIT1(a) \
+		case (a+1): SETBIT;
+#define SETBIT2(a) \
+		SETBIT1(a+1); \
+		SETBIT1(a);
+#define SETBIT4(a) \
+		SETBIT2(a+2); \
+		SETBIT2(a);
+#define SETBIT8(a) \
+		SETBIT4(a+4); \
+		SETBIT4(a);
+#define SETBIT16(a) \
+		SETBIT8(a+8); \
+		SETBIT8(a);
+#define SETBIT32(a) \
+		SETBIT16(a+16); \
+		SETBIT16(a);
+
+		int j = 0;
+		FBV_UINT jm = 1;
+		const FBV_UINT *lr = l.racks();
+		while ( j < n )
+		{
+			if ( jm == 0 )
+			{
+				jm = 1;
+				++lr;
+			}
+			
+			int dj = n - j;
+			switch ( n - j )
+			{
+				default: dj = 32;
+				SETBIT32(0);
+			}
+			j += dj;
+		}
+		
+		return;
+	}
+
+};
+
+
+#endif