Simplify bit operations and support integral points

6 files changed, 115 insertions, 285 deletions

diff --git a/Hilbert/Algorithm.hpp b/Hilbert/Algorithm.hpp
index c79f700..85640fa 100644
--- a/Hilbert/Algorithm.hpp
+++ b/Hilbert/Algorithm.hpp
@@ -237,7 +237,7 @@ namespace Hilbert
 		d = D0;
 		l.reset();
 		for ( j = 0; j < n; j++ )
-			p[j].reset();
+			p[j] = 0;
 
 		ho = m*n;
 		
diff --git a/Hilbert/GetBits.hpp b/Hilbert/GetBits.hpp
index 461690c..9f7feb6 100644
--- a/Hilbert/GetBits.hpp
+++ b/Hilbert/GetBits.hpp
@@ -1,4 +1,4 @@
-/*
+/* 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
@@ -19,74 +19,21 @@
 #ifndef _GETBITS_HPP_
 #define _GETBITS_HPP_
 
+#include <Hilbert/Operations.hpp>
 
-#include <Hilbert/BigBitVec.hpp>
+namespace Hilbert {
 
-
-namespace Hilbert
+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
 {
-	template <class H,class I>
-	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.set(j,h.test(i+j));
-		return;
-	}
-
-	
-	// <CBigBitVec,CBigBitVec>
-	// #TODO
-
-	// <CBigBitVec,CFixBitVec>
-	template<>
-	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<>
-	inline
-	void
-	getBits(
-		const CFixBitVec &h,
-		int n,
-		int i,
-		CFixBitVec &w
-		)
-	{
-		w = h >> i;
-		w.truncate(n);
-		return;
+	for (int j = 0; j < n; j++) {
+		setBit(w, j, testBit(h, i + j));
 	}
+}
 
-};
-
+} // namespace Hilbert
 
 #endif
diff --git a/Hilbert/GetLocation.hpp b/Hilbert/GetLocation.hpp
index e66dc04..7fc2bca 100644
--- a/Hilbert/GetLocation.hpp
+++ b/Hilbert/GetLocation.hpp
@@ -1,4 +1,4 @@
-/*
+/* 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
@@ -19,95 +19,19 @@
 #ifndef _GETLOCATION_HPP_
 #define _GETLOCATION_HPP_
 
+#include <Hilbert/Operations.hpp>
 
-#include <Hilbert/BigBitVec.hpp>
+namespace Hilbert {
 
-
-namespace Hilbert
+template <class P, class I>
+inline void
+getLocation(const P* const p, const int n, const int i, I& l)
 {
-	
-	template<class P,class I>
-	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>
-	inline
-	void
-	getLocation(
-		const P *p,
-		int n,
-		int i,
-		I	&l
-		)
-	{
-		/*int j;
-		for ( j = n-1; j >= 0; --j )
-			l.set(j,p[j].test(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;
+	for (int j = 0; j < n; ++j) {
+		setBit(l, j, testBit(p[j], i));
 	}
-	
-};
+}
 
+} // namespace Hilbert
 
 #endif
diff --git a/Hilbert/Operations.hpp b/Hilbert/Operations.hpp
new file mode 100644
index 0000000..1c76dc2
--- /dev/null
+++ b/Hilbert/Operations.hpp
@@ -0,0 +1,70 @@
+/*
+ * Copyright (C) 2018 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 _OPERATIONS_HPP_
+#define _OPERATIONS_HPP_
+
+#include <cassert>
+#include <climits>
+#include <cstddef>
+#include <type_traits>
+
+/// 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);
+}
+
+#endif
diff --git a/Hilbert/SetBits.hpp b/Hilbert/SetBits.hpp
index 5f311f1..a040b8a 100644
--- a/Hilbert/SetBits.hpp
+++ b/Hilbert/SetBits.hpp
@@ -1,4 +1,4 @@
-/*
+/* 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
@@ -19,74 +19,21 @@
 #ifndef _SETBITS_HPP_
 #define _SETBITS_HPP_
 
+#include <Hilbert/Operations.hpp>
 
-#include <Hilbert/BigBitVec.hpp>
+namespace Hilbert {
 
-
-namespace Hilbert
+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
 {
-
-	template<class H,class I>
-	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.set(i+j,w.test(j));
+	for (int j = 0; j < n; j++) {
+		setBit(h, i + j, testBit(w, j));
 	}
+}
 
-
-	// <CBigBitVec,CBigBitVec>
-	// #TODO
-
-
-	// <CBigBitVec,CFixBitVec>
-	template<>
-	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<>
-	inline
-	void
-	setBits(
-		CFixBitVec &h,
-		int n,
-		int i,
-		const CFixBitVec &w
-		)
-	{
-		h |= w << i;
-		return;
-	}
-
-};
-
+} // namespace Hilbert
 
 #endif
diff --git a/Hilbert/SetLocation.hpp b/Hilbert/SetLocation.hpp
index 66e6459..4a2ba3d 100644
--- a/Hilbert/SetLocation.hpp
+++ b/Hilbert/SetLocation.hpp
@@ -1,4 +1,4 @@
-/*
+/* 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
@@ -19,77 +19,19 @@
 #ifndef _SETLOCATION_HPP_
 #define _SETLOCATION_HPP_
 
+#include <Hilbert/Operations.hpp>
 
-#include <Hilbert/BigBitVec.hpp>
+namespace Hilbert {
 
-
-namespace Hilbert
+template <class P, class I>
+inline void
+setLocation(P* const p, const int n, const int i, const I& l)
 {
-	template <class P,class I>
-	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].set(i,l.test(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;
+	for (int j = 0; j < n; j++) {
+		setBit(p[j], i, testBit(l, j));
 	}
+}
 
-};
-
+} // namespace Hilbert
 
 #endif