1 /** 2 * @file hilbert.h 3 * 4 * @section LICENSE 5 * 6 * The MIT License 7 * 8 * @copyright Copyright (c) 2020-2021 TileDB, Inc. 9 * 10 * Permission is hereby granted, free of charge, to any person obtaining a copy 11 * of this software and associated documentation files (the "Software"), to deal 12 * in the Software without restriction, including without limitation the rights 13 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 14 * copies of the Software, and to permit persons to whom the Software is 15 * furnished to do so, subject to the following conditions: 16 * 17 * The above copyright notice and this permission notice shall be included in 18 * all copies or substantial portions of the Software. 19 * 20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 21 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 23 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 25 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 26 * THE SOFTWARE. 27 * 28 * @section DESCRIPTION 29 * 30 * This file defines class Hilbert. This class utilizes two functions from 31 * John Skilling, "Programming the Hilbert Curve". In AIP, 2004. 32 */ 33 34 #ifndef TILEDB_HILBERT_H 35 #define TILEDB_HILBERT_H 36 37 #include <sys/types.h> 38 #include <cmath> 39 #include <cstdlib> 40 #include <iostream> 41 42 namespace tiledb { 43 namespace sm { 44 45 /* ********************************* */ 46 /* CONSTANTS */ 47 /* ********************************* */ 48 49 /** 50 * The Hilbert curve fills a multi-dimensional space in a particular manner 51 * with a 1D line. The typical operations of this class is converting a 52 * multi-dimensional tuple of coordinates into a 1D Hilbert value, and 53 * vice versa. 54 * 55 * For the 2D case, the Hilbert curve looks as follows: 56 * 57 \verbatim 58 | 59 15 | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ 60 | | | | | | | | | | | | | | | | | 61 | @ @---@ @ @ @---@ @ @ @---@ @ @ @---@ @ 62 | | | | | | | | | 63 | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ 64 | | | | | | | | | 65 | @---@ @---@---@---@ @---@ @---@ @---@---@---@ @---@ 66 | | | | | 67 | @ @---@---@ @---@---@ @ @ @---@---@ @---@---@ @ 68 | | | | | | | | | | | | | 69 Dim[1]| @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ 70 | | | | | 71 | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ 72 | | | | | | | | | | | | | 73 | @ @---@---@ @---@---@ @---@ @---@---@ @---@---@ @ 74 | | | 75 | @---@ @---@---@ @---@---@ @---@---@ @---@---@ @---@ 76 | | | | | | | | | | | 77 | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@ 78 | | | | | | | 79 | @ @---@ @ @---@ @---@ @---@ @---@ @ @---@ @ 80 | | | | | | | | | | | | | | | 81 | @---@ @---@ @ @---@---@ @---@---@ @ @---@ @---@ 82 | | | 83 3 | 5---6 9---@ @ @---@---@ @---@---@ @ @---@ @---@ 84 | | | | | | | | | | | | | | | 85 2 | 4 7---8 @ 30---@ @---@ @---@ @---@ @ @---@ @ 86 | | | | | | | 87 1 | 3---2 @---@ @---@ @---@ @---@ @---@ @---@ @---@ 88 | | | | | | | | | | | 89 0 | 0---1 @---@---@ @--20---@ @---@---@ @---@---@ @--255 90 | 91 ------------------------------------------------------------------- 92 0 1 2 3 Dim[0] 15 93 94 \endverbatim 95 * 96 * The Hilbert value of (2,3) is 9, whereas the coordinates corresponding to 97 * Hilbert value 2 are (1,1). 98 * 99 * The class utilizes two functions from John Skilling's work published in: 100 * John Skilling, "Programming the Hilbert Curve". In AIP, 2004 101 */ 102 class Hilbert { 103 public: 104 /* ********************************* */ 105 /* CONSTANTS */ 106 /* ********************************* */ 107 108 /** 109 * Maximum number of dimensions for defining the Hilbert curve. Although the 110 * Hilbert curve can be defined over arbitrary dimensionality, we limit the 111 * number of dimensions because they affect the number of bits used to 112 * represent a Hilbert value; in our class, a Hilbert value is a int64_t 113 * number and, thus, it cannot exceed 64 bits. 114 */ 115 static const int HC_MAX_DIM = 16; 116 117 /* ********************************* */ 118 /* CONSTRUCTORS & DESTRUCTORS */ 119 /* ********************************* */ 120 121 /** 122 * Constructor. 123 * 124 * @param bits Number of bits used for coordinate values across each dimension 125 * @param dim_num Number of dimensions 126 */ Hilbert(int bits,int dim_num)127 Hilbert(int bits, int dim_num) 128 : bits_(bits) 129 , dim_num_(dim_num) { 130 assert(dim_num >= 0 && dim_num < HC_MAX_DIM); 131 assert(bits * dim_num <= int(sizeof(uint64_t) * 8 - 1)); 132 } 133 134 /** 135 * Constructor. The number of bits will be calculated based 136 * on the number of dimensions. 137 * 138 * @param dim_num Number of dimensions 139 */ Hilbert(int dim_num)140 Hilbert(int dim_num) 141 : dim_num_(dim_num) { 142 assert(dim_num >= 0 && dim_num < HC_MAX_DIM); 143 bits_ = (sizeof(uint64_t) * 8 - 1) / dim_num_; 144 } 145 146 /** Destructor. */ 147 ~Hilbert() = default; 148 149 /* ********************************* */ 150 /* API */ 151 /* ********************************* */ 152 153 /** Returns the number of bits. */ bits()154 int bits() const { 155 return bits_; 156 } 157 158 /** Returns the number of dimensions. */ dim_num()159 int dim_num() const { 160 return dim_num_; 161 } 162 163 /** 164 * Converts a set of coordinates to a Hilbert value. 165 * 166 * @param coords The coordinates to be converted. 167 * @return The output Hilbert value. 168 */ coords_to_hilbert(uint64_t * coords)169 uint64_t coords_to_hilbert(uint64_t* coords) { 170 assert(coords != nullptr); 171 uint64_t ret = 0; 172 173 // Convert coords to the transpose form of the hilbert value 174 axes_to_transpose(coords, bits_, dim_num_); 175 176 // Convert the hilbert transpose form into an int64_t hilbert value 177 uint64_t c = 1; // This is a bit shifted from right to left over coords[i] 178 uint64_t h = 1; // This is a bit shifted from right to left over hilbert 179 for (int j = 0; j < bits_; ++j, c <<= 1) { 180 for (int i = dim_num_ - 1; i >= 0; --i, h <<= 1) { 181 if (coords[i] & c) 182 ret |= h; 183 } 184 } 185 186 return ret; 187 } 188 189 /** 190 * Converts a Hilbert value into a set of coordinates. 191 * 192 * @param hilbert The Hilbert value to be converted. 193 * @param coords The output coordinates. 194 * @return void 195 */ hilbert_to_coords(uint64_t hilbert,uint64_t * coords)196 void hilbert_to_coords(uint64_t hilbert, uint64_t* coords) { 197 assert(coords != nullptr); 198 199 // Initialization 200 for (int i = 0; i < dim_num_; ++i) 201 coords[i] = 0; 202 203 // Convert the int64_t hilbert value to its transpose form 204 uint64_t c = 1; // This is a bit shifted from right to left over coords[i] 205 uint64_t h = 1; // This is a bit shifted from right to left over hilbert 206 for (int j = 0; j < bits_; ++j, c <<= 1) { 207 for (int i = dim_num_ - 1; i >= 0; --i, h <<= 1) { 208 if (hilbert & h) 209 coords[i] |= c; 210 } 211 } 212 213 // Convert coords to the transpose form of the hilbert value 214 transpose_to_axes(coords, bits_, dim_num_); 215 } 216 217 private: 218 /* ********************************* */ 219 /* PRIVATE ATTRIBUTES */ 220 /* ********************************* */ 221 222 /** Number of bits for representing a coordinate per dimension. */ 223 int bits_; 224 /** Number of dimensions. */ 225 int dim_num_; 226 227 /* ********************************* */ 228 /* PRIVATE METHODS */ 229 /* ********************************* */ 230 231 /** 232 * From John Skilling's work. It converts the input coordinates 233 * to what is called the transpose of the Hilbert value. This is done 234 * in place. 235 * 236 * **Example** 237 * 238 * Let b=5 and n=3. Let the 15-bit Hilbert value of the input coordinates 239 * be A B C D E a b c d e 1 2 3 4 5. The function places this number into 240 * parameter X as follows: 241 * 242 \verbatim 243 X[0] = A D b e 3 X[1]| 244 X[1] = B E c 1 4 <-------> | /X[2] 245 X[2] = C a d 2 5 axes | / 246 high low |/______ 247 248 X[0] 249 \endverbatim 250 * 251 * The X value after the function terminates is called the transpose form 252 * of the Hilbert value. 253 * 254 * @param X Input coordinates, and output transpose (the conversion is 255 * done in place). 256 * @param b Number of bits for representing a coordinate per dimension 257 * (it should be equal to bits_). 258 * @param n Number of dimensions (it should be equal to dim_num_). 259 * @return void 260 */ axes_to_transpose(uint64_t * X,int b,int n)261 void axes_to_transpose(uint64_t* X, int b, int n) { 262 uint64_t P, Q, t; 263 int i; 264 265 // Inverse undo 266 for (Q = (uint64_t)1 << (b - 1); Q > 1; Q >>= 1) { 267 P = Q - 1; 268 if (X[0] & Q) // invert 269 X[0] ^= P; 270 for (i = 1; i < n; i++) 271 if (X[i] & Q) // invert 272 X[0] ^= P; 273 else { // exchange 274 t = (X[0] ^ X[i]) & P; 275 X[0] ^= t; 276 X[i] ^= t; 277 } 278 } 279 280 // Gray encode (inverse of decode) 281 for (i = 1; i < n; i++) 282 X[i] ^= X[i - 1]; 283 t = X[n - 1]; 284 for (i = 1; i < b; i <<= 1) 285 X[n - 1] ^= X[n - 1] >> i; 286 t ^= X[n - 1]; 287 for (i = n - 2; i >= 0; i--) 288 X[i] ^= t; 289 } 290 291 /** 292 * From John Skilling's work. It converts the transpose of a 293 * Hilbert value into the corresponding coordinates. This is done in place. 294 * 295 * **Example** 296 * 297 * Let b=5 and n=3. Let the 15-bit Hilbert value of the output coordinates 298 * be A B C D E a b c d e 1 2 3 4 5. The function takes as input the tranpose 299 * form of the Hilbert value, which is stored in X as follows: 300 * 301 \verbatim 302 X[0] = A D b e 3 X[1]| 303 X[1] = B E c 1 4 <-------> | /X[2] 304 X[2] = C a d 2 5 axes | / 305 high low |/______ 306 307 X[0] 308 \endverbatim 309 * 310 * The X value after the function terminates will contain the coordinates 311 * corresponding to the Hilbert value. 312 * 313 * @param X Input transpose, and output coordinates (the conversion is 314 * done in place). 315 * @param b Number of bits for representing a coordinate per dimension 316 * (it should be equal to bits_). 317 * @param n Number of dimensions (it should be equal to dim_num_). 318 * @return void 319 */ transpose_to_axes(uint64_t * X,int b,int n)320 void transpose_to_axes(uint64_t* X, int b, int n) { 321 uint64_t M, P, Q, t; 322 int i; 323 324 // Gray decode by H ^ (H/2) 325 t = X[n - 1] >> 1; 326 for (i = n - 1; i; i--) 327 X[i] ^= X[i - 1]; 328 X[0] ^= t; 329 330 // Undo excess work 331 M = 2 << (b - 1); 332 for (Q = 2; Q != M; Q <<= 1) { 333 P = Q - 1; 334 for (i = n - 1; i; i--) 335 if (X[i] & Q) // invert 336 X[0] ^= P; 337 else { // exchange 338 t = (X[0] ^ X[i]) & P; 339 X[0] ^= t; 340 X[i] ^= t; 341 } 342 if (X[0] & Q) // invert 343 X[0] ^= P; 344 } 345 } 346 }; 347 348 } // namespace sm 349 } // end namespace tiledb 350 351 #endif // TILEDB_HILBERT_H 352