diff options
| -rw-r--r-- | Hilbert/Algorithm.hpp | 482 | ||||
| -rw-r--r-- | Hilbert/BigBitVec.hpp | 928 | ||||
| -rw-r--r-- | Hilbert/Common.hpp | 26 | ||||
| -rw-r--r-- | Hilbert/FixBitVec.hpp | 538 | ||||
| -rw-r--r-- | Hilbert/GetBits.hpp | 93 | ||||
| -rw-r--r-- | Hilbert/GetLocation.hpp | 114 | ||||
| -rw-r--r-- | Hilbert/GrayCodeRank.hpp | 246 | ||||
| -rw-r--r-- | Hilbert/Hilbert.hpp | 142 | ||||
| -rw-r--r-- | Hilbert/SetBits.hpp | 93 | ||||
| -rw-r--r-- | Hilbert/SetLocation.hpp | 96 |
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 |