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