1 |
greg |
3.2 |
#ifndef lint |
2 |
|
|
static const char RCSid[] = "$Id$"; |
3 |
|
|
#endif |
4 |
greg |
3.1 |
/* See LICENSE below for information on rights to use, modify and distribute |
5 |
|
|
this code. */ |
6 |
|
|
|
7 |
|
|
/* |
8 |
|
|
* hilbert.c - Computes Hilbert space-filling curve coordinates, without |
9 |
|
|
* recursion, from integer index, and vice versa, and other Hilbert-related |
10 |
|
|
* calculations. Also known as Pi-order or Peano scan. |
11 |
|
|
* |
12 |
|
|
* Author: Doug Moore |
13 |
|
|
* Dept. of Computational and Applied Math |
14 |
|
|
* Rice University |
15 |
|
|
* http://www.caam.rice.edu/~dougm |
16 |
|
|
* Date: Sun Feb 20 2000 |
17 |
|
|
* Copyright (c) 1998-2000, Rice University |
18 |
|
|
* |
19 |
|
|
* Acknowledgement: |
20 |
|
|
* This implementation is based on the work of A. R. Butz ("Alternative |
21 |
|
|
* Algorithm for Hilbert's Space-Filling Curve", IEEE Trans. Comp., April, |
22 |
|
|
* 1971, pp 424-426) and its interpretation by Spencer W. Thomas, University |
23 |
|
|
* of Michigan (http://www-personal.umich.edu/~spencer/Home.html) in his widely |
24 |
|
|
* available C software. While the implementation here differs considerably |
25 |
|
|
* from his, the first two interfaces and the style of some comments are very |
26 |
|
|
* much derived from his work. */ |
27 |
|
|
|
28 |
|
|
|
29 |
|
|
#include "hilbert.h" |
30 |
|
|
|
31 |
|
|
/* implementation of the hilbert functions */ |
32 |
|
|
|
33 |
|
|
#define adjust_rotation(rotation,nDims,bits) \ |
34 |
|
|
do { \ |
35 |
|
|
/* rotation = (rotation + 1 + ffs(bits)) % nDims; */ \ |
36 |
|
|
bits &= -bits & nd1Ones; \ |
37 |
|
|
while (bits) \ |
38 |
|
|
bits >>= 1, ++rotation; \ |
39 |
|
|
if ( ++rotation >= nDims ) \ |
40 |
|
|
rotation -= nDims; \ |
41 |
|
|
} while (0) |
42 |
|
|
|
43 |
|
|
#define ones(T,k) ((((T)2) << (k-1)) - 1) |
44 |
|
|
|
45 |
|
|
#define rdbit(w,k) (((w) >> (k)) & 1) |
46 |
|
|
|
47 |
|
|
#define rotateRight(arg, nRots, nDims) \ |
48 |
|
|
((((arg) >> (nRots)) | ((arg) << ((nDims)-(nRots)))) & ones(bitmask_t,nDims)) |
49 |
|
|
|
50 |
|
|
#define rotateLeft(arg, nRots, nDims) \ |
51 |
|
|
((((arg) << (nRots)) | ((arg) >> ((nDims)-(nRots)))) & ones(bitmask_t,nDims)) |
52 |
|
|
|
53 |
|
|
#define DLOGB_BIT_TRANSPOSE |
54 |
|
|
static bitmask_t |
55 |
|
|
bitTranspose(unsigned nDims, unsigned nBits, bitmask_t inCoords) |
56 |
|
|
#if defined(DLOGB_BIT_TRANSPOSE) |
57 |
|
|
{ |
58 |
|
|
unsigned const nDims1 = nDims-1; |
59 |
|
|
unsigned inB = nBits; |
60 |
|
|
unsigned utB; |
61 |
|
|
bitmask_t inFieldEnds = 1; |
62 |
|
|
bitmask_t inMask = ones(bitmask_t,inB); |
63 |
|
|
bitmask_t coords = 0; |
64 |
|
|
|
65 |
|
|
while ((utB = inB >> 1)) |
66 |
|
|
{ |
67 |
|
|
unsigned const shiftAmt = nDims1 * utB; |
68 |
|
|
bitmask_t const utFieldEnds = |
69 |
|
|
inFieldEnds | (inFieldEnds << (shiftAmt+utB)); |
70 |
|
|
bitmask_t const utMask = |
71 |
|
|
(utFieldEnds << utB) - utFieldEnds; |
72 |
|
|
bitmask_t utCoords = 0; |
73 |
|
|
unsigned d; |
74 |
|
|
if (inB & 1) |
75 |
|
|
{ |
76 |
|
|
bitmask_t const inFieldStarts = inFieldEnds << (inB-1); |
77 |
|
|
unsigned oddShift = 2*shiftAmt; |
78 |
|
|
for (d = 0; d < nDims; ++d) |
79 |
|
|
{ |
80 |
|
|
bitmask_t in = inCoords & inMask; |
81 |
|
|
inCoords >>= inB; |
82 |
|
|
coords |= (in & inFieldStarts) << oddShift++; |
83 |
|
|
in &= ~inFieldStarts; |
84 |
|
|
in = (in | (in << shiftAmt)) & utMask; |
85 |
|
|
utCoords |= in << (d*utB); |
86 |
|
|
} |
87 |
|
|
} |
88 |
|
|
else |
89 |
|
|
{ |
90 |
|
|
for (d = 0; d < nDims; ++d) |
91 |
|
|
{ |
92 |
|
|
bitmask_t in = inCoords & inMask; |
93 |
|
|
inCoords >>= inB; |
94 |
|
|
in = (in | (in << shiftAmt)) & utMask; |
95 |
|
|
utCoords |= in << (d*utB); |
96 |
|
|
} |
97 |
|
|
} |
98 |
|
|
inCoords = utCoords; |
99 |
|
|
inB = utB; |
100 |
|
|
inFieldEnds = utFieldEnds; |
101 |
|
|
inMask = utMask; |
102 |
|
|
} |
103 |
|
|
coords |= inCoords; |
104 |
|
|
return coords; |
105 |
|
|
} |
106 |
|
|
#else |
107 |
|
|
{ |
108 |
|
|
bitmask_t coords = 0; |
109 |
|
|
unsigned d; |
110 |
|
|
for (d = 0; d < nDims; ++d) |
111 |
|
|
{ |
112 |
|
|
unsigned b; |
113 |
|
|
bitmask_t in = inCoords & ones(bitmask_t,nBits); |
114 |
|
|
bitmask_t out = 0; |
115 |
|
|
inCoords >>= nBits; |
116 |
|
|
for (b = nBits; b--;) |
117 |
|
|
{ |
118 |
|
|
out <<= nDims; |
119 |
|
|
out |= rdbit(in, b); |
120 |
|
|
} |
121 |
|
|
coords |= out << d; |
122 |
|
|
} |
123 |
|
|
return coords; |
124 |
|
|
} |
125 |
|
|
#endif |
126 |
|
|
|
127 |
|
|
/***************************************************************** |
128 |
|
|
* hilbert_i2c |
129 |
|
|
* |
130 |
|
|
* Convert an index into a Hilbert curve to a set of coordinates. |
131 |
|
|
* Inputs: |
132 |
|
|
* nDims: Number of coordinate axes. |
133 |
|
|
* nBits: Number of bits per axis. |
134 |
|
|
* index: The index, contains nDims*nBits bits |
135 |
|
|
* (so nDims*nBits must be <= 8*sizeof(bitmask_t)). |
136 |
|
|
* Outputs: |
137 |
|
|
* coord: The list of nDims coordinates, each with nBits bits. |
138 |
|
|
* Assumptions: |
139 |
|
|
* nDims*nBits <= (sizeof index) * (bits_per_byte) |
140 |
|
|
*/ |
141 |
|
|
void |
142 |
|
|
hilbert_i2c(unsigned nDims, unsigned nBits, bitmask_t index, bitmask_t coord[]) |
143 |
|
|
{ |
144 |
|
|
if (nDims > 1) |
145 |
|
|
{ |
146 |
|
|
bitmask_t coords; |
147 |
|
|
halfmask_t const nbOnes = ones(halfmask_t,nBits); |
148 |
|
|
unsigned d; |
149 |
|
|
|
150 |
|
|
if (nBits > 1) |
151 |
|
|
{ |
152 |
|
|
unsigned const nDimsBits = nDims*nBits; |
153 |
|
|
halfmask_t const ndOnes = ones(halfmask_t,nDims); |
154 |
|
|
halfmask_t const nd1Ones= ndOnes >> 1; /* for adjust_rotation */ |
155 |
|
|
unsigned b = nDimsBits; |
156 |
|
|
unsigned rotation = 0; |
157 |
|
|
halfmask_t flipBit = 0; |
158 |
|
|
bitmask_t const nthbits = ones(bitmask_t,nDimsBits) / ndOnes; |
159 |
|
|
index ^= (index ^ nthbits) >> 1; |
160 |
|
|
coords = 0; |
161 |
|
|
do |
162 |
|
|
{ |
163 |
|
|
halfmask_t bits = (index >> (b-=nDims)) & ndOnes; |
164 |
|
|
coords <<= nDims; |
165 |
|
|
coords |= rotateLeft(bits, rotation, nDims) ^ flipBit; |
166 |
|
|
flipBit = (halfmask_t)1 << rotation; |
167 |
|
|
adjust_rotation(rotation,nDims,bits); |
168 |
|
|
} while (b); |
169 |
|
|
for (b = nDims; b < nDimsBits; b *= 2) |
170 |
|
|
coords ^= coords >> b; |
171 |
|
|
coords = bitTranspose(nBits, nDims, coords); |
172 |
|
|
} |
173 |
|
|
else |
174 |
|
|
coords = index ^ (index >> 1); |
175 |
|
|
|
176 |
|
|
for (d = 0; d < nDims; ++d) |
177 |
|
|
{ |
178 |
|
|
coord[d] = coords & nbOnes; |
179 |
|
|
coords >>= nBits; |
180 |
|
|
} |
181 |
|
|
} |
182 |
|
|
else |
183 |
|
|
coord[0] = index; |
184 |
|
|
} |
185 |
|
|
|
186 |
|
|
/***************************************************************** |
187 |
|
|
* hilbert_c2i |
188 |
|
|
* |
189 |
|
|
* Convert coordinates of a point on a Hilbert curve to its index. |
190 |
|
|
* Inputs: |
191 |
|
|
* nDims: Number of coordinates. |
192 |
|
|
* nBits: Number of bits/coordinate. |
193 |
|
|
* coord: Array of n nBits-bit coordinates. |
194 |
|
|
* Outputs: |
195 |
|
|
* index: Output index value. nDims*nBits bits. |
196 |
|
|
* Assumptions: |
197 |
|
|
* nDims*nBits <= (sizeof bitmask_t) * (bits_per_byte) |
198 |
|
|
*/ |
199 |
|
|
bitmask_t |
200 |
|
|
hilbert_c2i(unsigned nDims, unsigned nBits, bitmask_t const coord[]) |
201 |
|
|
{ |
202 |
|
|
if (nDims > 1) |
203 |
|
|
{ |
204 |
|
|
unsigned const nDimsBits = nDims*nBits; |
205 |
|
|
bitmask_t index; |
206 |
|
|
unsigned d; |
207 |
|
|
bitmask_t coords = 0; |
208 |
|
|
for (d = nDims; d--; ) |
209 |
|
|
{ |
210 |
|
|
coords <<= nBits; |
211 |
|
|
coords |= coord[d]; |
212 |
|
|
} |
213 |
|
|
|
214 |
|
|
if (nBits > 1) |
215 |
|
|
{ |
216 |
|
|
halfmask_t const ndOnes = ones(halfmask_t,nDims); |
217 |
|
|
halfmask_t const nd1Ones= ndOnes >> 1; /* for adjust_rotation */ |
218 |
|
|
unsigned b = nDimsBits; |
219 |
|
|
unsigned rotation = 0; |
220 |
|
|
halfmask_t flipBit = 0; |
221 |
|
|
bitmask_t const nthbits = ones(bitmask_t,nDimsBits) / ndOnes; |
222 |
|
|
coords = bitTranspose(nDims, nBits, coords); |
223 |
|
|
coords ^= coords >> nDims; |
224 |
|
|
index = 0; |
225 |
|
|
do |
226 |
|
|
{ |
227 |
|
|
halfmask_t bits = (coords >> (b-=nDims)) & ndOnes; |
228 |
|
|
bits = rotateRight(flipBit ^ bits, rotation, nDims); |
229 |
|
|
index <<= nDims; |
230 |
|
|
index |= bits; |
231 |
|
|
flipBit = (halfmask_t)1 << rotation; |
232 |
|
|
adjust_rotation(rotation,nDims,bits); |
233 |
|
|
} while (b); |
234 |
|
|
index ^= nthbits >> 1; |
235 |
|
|
} |
236 |
|
|
else |
237 |
|
|
index = coords; |
238 |
|
|
for (d = 1; d < nDimsBits; d *= 2) |
239 |
|
|
index ^= index >> d; |
240 |
|
|
return index; |
241 |
|
|
} |
242 |
|
|
else |
243 |
|
|
return coord[0]; |
244 |
|
|
} |
245 |
|
|
|
246 |
|
|
/***************************************************************** |
247 |
|
|
* Readers and writers of bits |
248 |
|
|
*/ |
249 |
|
|
|
250 |
|
|
typedef bitmask_t (*BitReader) (unsigned nDims, unsigned nBytes, |
251 |
|
|
char const* c, unsigned y); |
252 |
|
|
typedef void (*BitWriter) (unsigned d, unsigned nBytes, |
253 |
|
|
char* c, unsigned y, int fold); |
254 |
|
|
|
255 |
|
|
|
256 |
|
|
#if defined(sparc) |
257 |
|
|
#define __BIG_ENDIAN__ |
258 |
|
|
#endif |
259 |
|
|
|
260 |
|
|
#if defined(__BIG_ENDIAN__) |
261 |
|
|
#define whichByte(nBytes,y) (nBytes-1-y/8) |
262 |
|
|
#define setBytes(dst,pos,nBytes,val) \ |
263 |
|
|
memset(&dst[pos+1],val,nBytes-pos-1) |
264 |
|
|
#else |
265 |
|
|
#define whichByte(nBytes,y) (y/8) |
266 |
|
|
#define setBytes(dst,pos,nBytes,val) \ |
267 |
|
|
memset(&dst[0],val,pos) |
268 |
|
|
#endif |
269 |
|
|
|
270 |
|
|
static bitmask_t |
271 |
|
|
getIntBits(unsigned nDims, unsigned nBytes, char const* c, unsigned y) |
272 |
|
|
{ |
273 |
|
|
unsigned const bit = y%8; |
274 |
|
|
unsigned const offs = whichByte(nBytes,y); |
275 |
|
|
unsigned d; |
276 |
|
|
bitmask_t bits = 0; |
277 |
|
|
c += offs; |
278 |
|
|
for (d = 0; d < nDims; ++d) |
279 |
|
|
{ |
280 |
|
|
bits |= rdbit(*c, bit) << d; |
281 |
|
|
c += nBytes; |
282 |
|
|
} |
283 |
|
|
return bits; |
284 |
|
|
} |
285 |
|
|
|
286 |
|
|
#include <string.h> |
287 |
|
|
static void |
288 |
|
|
propogateIntBits(unsigned d, unsigned nBytes, |
289 |
|
|
char* c, unsigned y, int fold) |
290 |
|
|
{ |
291 |
|
|
unsigned const byteId = whichByte(nBytes,y); |
292 |
|
|
unsigned const b = y%8; |
293 |
|
|
char const bthbit = 1 << b; |
294 |
|
|
char* const target = &c[d*nBytes]; |
295 |
|
|
target[byteId] ^= bthbit; |
296 |
|
|
if (!fold) |
297 |
|
|
{ |
298 |
|
|
char notbit = ((target[byteId] >> b) & 1) - 1; |
299 |
|
|
if (notbit) |
300 |
|
|
target[byteId] |= bthbit-1; |
301 |
|
|
else |
302 |
|
|
target[byteId] &= -bthbit; |
303 |
|
|
setBytes(target,byteId,nBytes,notbit); |
304 |
|
|
} |
305 |
|
|
} |
306 |
|
|
|
307 |
|
|
/* An IEEE double is treated as a 2100 bit number. In particular, 0 is treated |
308 |
|
|
as a 1 followed by 2099 zeroes, and negative 0 as a 0 followed by 2099 ones. |
309 |
|
|
Only 53 bits differ between a number and a zero of the same sign, with the |
310 |
|
|
position of the 53 determined by the exponent, and the values of the 53 by |
311 |
|
|
the significand (with implicit leading 1 bit). Although IEEE 754 uses the |
312 |
|
|
maximum exponent for NaN's and infinities, this implementation ignores that |
313 |
|
|
decision, so that infinities and NaN's are treated as very large numbers. |
314 |
|
|
Note that we do not explicitly construct a 2100 bit bitmask in the IEEE |
315 |
|
|
routines below. */ |
316 |
|
|
|
317 |
|
|
enum { IEEEexpBits = 11 }; |
318 |
|
|
enum { IEEEsigBits = 52 }; |
319 |
|
|
enum { IEEErepBits = (1 << IEEEexpBits) + IEEEsigBits }; |
320 |
|
|
|
321 |
|
|
typedef union ieee754_double |
322 |
|
|
{ |
323 |
|
|
double d; |
324 |
|
|
|
325 |
|
|
/* This is the IEEE 754 double-precision format. */ |
326 |
|
|
struct |
327 |
|
|
{ |
328 |
|
|
#if defined(__BIG_ENDIAN__) |
329 |
|
|
unsigned int negative:1; |
330 |
|
|
unsigned int exponent:11; |
331 |
|
|
/* Together these comprise the mantissa. */ |
332 |
|
|
unsigned int mantissa0:20; |
333 |
|
|
unsigned int mantissa1:32; |
334 |
|
|
#else /* Big endian. */ |
335 |
|
|
/* Together these comprise the mantissa. */ |
336 |
|
|
unsigned int mantissa1:32; |
337 |
|
|
unsigned int mantissa0:20; |
338 |
|
|
unsigned int exponent:11; |
339 |
|
|
unsigned int negative:1; |
340 |
|
|
#endif /* Little endian. */ |
341 |
|
|
} ieee; |
342 |
|
|
} ieee754_double; |
343 |
|
|
|
344 |
|
|
static bitmask_t |
345 |
|
|
getIEEESignBits(unsigned nDims, double const* c) |
346 |
|
|
{ |
347 |
|
|
unsigned d; |
348 |
|
|
ieee754_double x; |
349 |
|
|
bitmask_t bits = 0; |
350 |
|
|
for (d = 0; d < nDims; ++d) |
351 |
|
|
{ |
352 |
|
|
x.d = c[d]; |
353 |
|
|
bits |= x.ieee.negative << d; |
354 |
|
|
} |
355 |
|
|
return bits; |
356 |
|
|
} |
357 |
|
|
|
358 |
|
|
static bitmask_t |
359 |
|
|
getIEEEBits(unsigned nDims, |
360 |
|
|
unsigned ignoreMe, /* ignored */ |
361 |
|
|
char const* cP, |
362 |
|
|
unsigned y) |
363 |
|
|
/* retrieve bits y of elements of double array c, where an expanded IEEE |
364 |
|
|
double has 2100 bits. */ |
365 |
|
|
{ |
366 |
|
|
unsigned d; |
367 |
|
|
double const* c = (double const*) cP; |
368 |
|
|
ieee754_double x; |
369 |
|
|
bitmask_t bits = 0; |
370 |
|
|
for (x.d = c[d=0]; d < nDims; x.d = c[++d]) |
371 |
|
|
{ |
372 |
|
|
bitmask_t bit = x.ieee.negative; |
373 |
|
|
unsigned normalized = (x.ieee.exponent != 0); |
374 |
|
|
unsigned diff = y - (x.ieee.exponent - normalized); |
375 |
|
|
if (diff <= 52) |
376 |
|
|
bit ^= 1 & ((diff < 32)? x.ieee.mantissa1 >> diff: |
377 |
|
|
(diff < 52)? x.ieee.mantissa0 >> (diff - 32): |
378 |
|
|
/* else */ normalized); |
379 |
|
|
else |
380 |
|
|
bit ^= (y == IEEErepBits-1); |
381 |
|
|
|
382 |
|
|
bits |= bit << d; |
383 |
|
|
} |
384 |
|
|
return bits; |
385 |
|
|
} |
386 |
|
|
|
387 |
|
|
static void |
388 |
|
|
propogateIEEEBits(unsigned d, unsigned nBytes, |
389 |
|
|
char* cP, unsigned y, int fold) |
390 |
|
|
{ |
391 |
|
|
ieee754_double* x = d + (ieee754_double*) cP; |
392 |
|
|
unsigned normalized = (x->ieee.exponent != 0); |
393 |
|
|
unsigned diff = y - (x->ieee.exponent - normalized); |
394 |
|
|
if (diff < 32) |
395 |
|
|
{ |
396 |
|
|
unsigned b = 1 << diff; |
397 |
|
|
unsigned bit = x->ieee.mantissa1 & b; |
398 |
|
|
x->ieee.mantissa1 &= ~(b-1); |
399 |
|
|
x->ieee.mantissa1 |= b; |
400 |
|
|
if (bit) |
401 |
|
|
--x->ieee.mantissa1; |
402 |
|
|
} |
403 |
|
|
else if (diff < 52) |
404 |
|
|
{ |
405 |
|
|
unsigned b = 1 << (diff - 32); |
406 |
|
|
unsigned bit = x->ieee.mantissa0 & b; |
407 |
|
|
x->ieee.mantissa0 &= ~(b-1); |
408 |
|
|
x->ieee.mantissa0 |= b; |
409 |
|
|
if (bit) |
410 |
|
|
--x->ieee.mantissa0; |
411 |
|
|
x->ieee.mantissa1 = bit?-1: 0; |
412 |
|
|
} |
413 |
|
|
else if (diff == 52) /* "flip" the implicit 1 bit */ |
414 |
|
|
{ |
415 |
|
|
if (normalized) |
416 |
|
|
--x->ieee.exponent; |
417 |
|
|
else |
418 |
|
|
x->ieee.exponent = 1; |
419 |
|
|
x->ieee.mantissa0 = -normalized; |
420 |
|
|
x->ieee.mantissa1 = -normalized; |
421 |
|
|
} |
422 |
|
|
else if (diff < IEEErepBits) |
423 |
|
|
{ |
424 |
|
|
if (y == IEEErepBits-1) |
425 |
|
|
{ |
426 |
|
|
x->ieee.negative ^= 1; |
427 |
|
|
x->ieee.exponent = 0; |
428 |
|
|
} |
429 |
|
|
else |
430 |
|
|
x->ieee.exponent = y - 51; |
431 |
|
|
x->ieee.mantissa0 = 0; |
432 |
|
|
x->ieee.mantissa1 = 0; |
433 |
|
|
} |
434 |
|
|
} |
435 |
|
|
|
436 |
|
|
static unsigned |
437 |
|
|
getIEEEexptMax(unsigned nDims, double const* c) |
438 |
|
|
{ |
439 |
|
|
unsigned max = 0; |
440 |
|
|
unsigned d; |
441 |
|
|
for (d = 0; d < nDims; ++d) |
442 |
|
|
{ |
443 |
|
|
ieee754_double x; |
444 |
|
|
x.d = c[d]; |
445 |
|
|
if (max < x.ieee.exponent) |
446 |
|
|
max = x.ieee.exponent; |
447 |
|
|
} |
448 |
|
|
if (max) --max; |
449 |
|
|
return max; |
450 |
|
|
} |
451 |
|
|
|
452 |
|
|
static void |
453 |
|
|
getIEEEinitValues(double const* c1, |
454 |
|
|
unsigned y, |
455 |
|
|
unsigned nDims, |
456 |
|
|
unsigned* rotation, |
457 |
|
|
bitmask_t* bits, |
458 |
|
|
bitmask_t* index) |
459 |
|
|
{ |
460 |
|
|
bitmask_t const one = 1; |
461 |
|
|
unsigned d; |
462 |
|
|
bitmask_t signBits = getIEEESignBits(nDims, c1); |
463 |
|
|
unsigned signParity, leastZeroBit, strayBit; |
464 |
|
|
|
465 |
|
|
/* compute the odd/evenness of the number of sign bits */ |
466 |
|
|
{ |
467 |
|
|
bitmask_t signPar = signBits; |
468 |
|
|
for (d = 1; d < nDims; d *= 2) |
469 |
|
|
signPar ^= signPar >> d; |
470 |
|
|
signParity = signPar & 1; |
471 |
|
|
} |
472 |
|
|
|
473 |
|
|
/* find the position of the least-order 0 bit in among signBits and adjust it |
474 |
|
|
if necessary */ |
475 |
|
|
for (leastZeroBit = 0; leastZeroBit < nDims; ++leastZeroBit) |
476 |
|
|
if (rdbit(signBits, leastZeroBit) == 0) |
477 |
|
|
break; |
478 |
|
|
strayBit = 0; |
479 |
|
|
if (leastZeroBit == nDims-2) |
480 |
|
|
strayBit = 1; |
481 |
|
|
else if (leastZeroBit == nDims) |
482 |
|
|
leastZeroBit = nDims-1; |
483 |
|
|
|
484 |
|
|
if (y % 2 == 1) |
485 |
|
|
{ |
486 |
|
|
*rotation = (IEEErepBits - y + 1 + leastZeroBit) % nDims; |
487 |
|
|
if (y < IEEErepBits-1) |
488 |
|
|
{ |
489 |
|
|
*bits = signBits ^ (one << ((*rotation + strayBit) % nDims)); |
490 |
|
|
*index = signParity; |
491 |
|
|
} |
492 |
|
|
else /* y == IEEErepBits-1 */ |
493 |
|
|
{ |
494 |
|
|
*bits = signBits ^ (ones(bitmask_t,nDims) &~ 1); |
495 |
|
|
*index = signParity ^ (nDims&1); |
496 |
|
|
} |
497 |
|
|
} |
498 |
|
|
else /* y % 2 == 0 */ |
499 |
|
|
if (y < IEEErepBits) |
500 |
|
|
{ |
501 |
|
|
unsigned shift_amt = (IEEErepBits - y + leastZeroBit) % nDims; |
502 |
|
|
*rotation = (shift_amt + 2 + strayBit) % nDims; |
503 |
|
|
*bits = signBits ^ (one << shift_amt); |
504 |
|
|
*index = signParity ^ 1; |
505 |
|
|
} |
506 |
|
|
else /* y == IEEErepBits */ |
507 |
|
|
{ |
508 |
|
|
*rotation = 0; |
509 |
|
|
*bits = one << (nDims-1); |
510 |
|
|
*index = 1; |
511 |
|
|
} |
512 |
|
|
} |
513 |
|
|
|
514 |
|
|
/***************************************************************** |
515 |
|
|
* hilbert_cmp, hilbert_ieee_cmp |
516 |
|
|
* |
517 |
|
|
* Determine which of two points lies further along the Hilbert curve |
518 |
|
|
* Inputs: |
519 |
|
|
* nDims: Number of coordinates. |
520 |
|
|
* nBytes: Number of bytes of storage/coordinate (hilbert_cmp only) |
521 |
|
|
* nBits: Number of bits/coordinate. (hilbert_cmp only) |
522 |
|
|
* coord1: Array of nDims nBytes-byte coordinates (or doubles for ieee_cmp). |
523 |
|
|
* coord2: Array of nDims nBytes-byte coordinates (or doubles for ieee_cmp). |
524 |
|
|
* Return value: |
525 |
|
|
* -1, 0, or 1 according to whether |
526 |
|
|
coord1<coord2, coord1==coord2, coord1>coord2 |
527 |
|
|
* Assumptions: |
528 |
|
|
* nBits <= (sizeof bitmask_t) * (bits_per_byte) |
529 |
|
|
*/ |
530 |
|
|
|
531 |
|
|
static int |
532 |
|
|
hilbert_cmp_work(unsigned nDims, unsigned nBytes, unsigned nBits, |
533 |
|
|
unsigned max, unsigned y, |
534 |
|
|
char const* c1, char const* c2, |
535 |
|
|
unsigned rotation, |
536 |
|
|
bitmask_t bits, |
537 |
|
|
bitmask_t index, |
538 |
|
|
BitReader getBits) |
539 |
|
|
{ |
540 |
|
|
bitmask_t const one = 1; |
541 |
|
|
bitmask_t const nd1Ones = ones(bitmask_t,nDims) >> 1; /* used in adjust_rotation macro */ |
542 |
|
|
while (y-- > max) |
543 |
|
|
{ |
544 |
|
|
bitmask_t reflection = getBits(nDims, nBytes, c1, y); |
545 |
|
|
bitmask_t diff = reflection ^ getBits(nDims, nBytes, c2, y); |
546 |
|
|
bits ^= reflection; |
547 |
|
|
bits = rotateRight(bits, rotation, nDims); |
548 |
|
|
if (diff) |
549 |
|
|
{ |
550 |
|
|
unsigned d; |
551 |
|
|
diff = rotateRight(diff, rotation, nDims); |
552 |
|
|
for (d = 1; d < nDims; d *= 2) |
553 |
|
|
{ |
554 |
|
|
index ^= index >> d; |
555 |
|
|
bits ^= bits >> d; |
556 |
|
|
diff ^= diff >> d; |
557 |
|
|
} |
558 |
|
|
return (((index ^ y ^ nBits) & 1) == (bits < (bits^diff)))? -1: 1; |
559 |
|
|
} |
560 |
|
|
index ^= bits; |
561 |
|
|
reflection ^= one << rotation; |
562 |
|
|
adjust_rotation(rotation,nDims,bits); |
563 |
|
|
bits = reflection; |
564 |
|
|
} |
565 |
|
|
return 0; |
566 |
|
|
} |
567 |
|
|
|
568 |
|
|
int |
569 |
|
|
hilbert_cmp(unsigned nDims, unsigned nBytes, unsigned nBits, |
570 |
|
|
void const* c1, void const* c2) |
571 |
|
|
{ |
572 |
|
|
bitmask_t const one = 1; |
573 |
|
|
bitmask_t bits = one << (nDims-1); |
574 |
|
|
return hilbert_cmp_work(nDims, nBytes, nBits, 0, nBits, |
575 |
|
|
(char const*)c1, (char const*)c2, |
576 |
|
|
0, bits, bits, getIntBits); |
577 |
|
|
} |
578 |
|
|
|
579 |
|
|
int |
580 |
|
|
hilbert_ieee_cmp(unsigned nDims, double const* c1, double const* c2) |
581 |
|
|
{ |
582 |
|
|
unsigned rotation, max; |
583 |
|
|
bitmask_t bits, index; |
584 |
|
|
if (getIEEESignBits(nDims, c1) != getIEEESignBits(nDims, c2)) |
585 |
|
|
max = 2047; |
586 |
|
|
else |
587 |
|
|
{ |
588 |
|
|
unsigned max1 = getIEEEexptMax(nDims, c1); |
589 |
|
|
unsigned max2 = getIEEEexptMax(nDims, c2); |
590 |
|
|
max = (max1 > max2)? max1: max2; |
591 |
|
|
} |
592 |
|
|
|
593 |
|
|
getIEEEinitValues(c1, max+53, nDims, &rotation, &bits, &index); |
594 |
|
|
return hilbert_cmp_work(nDims, 8, 64, max, max+53, |
595 |
|
|
(char const*)c1, (char const*)c2, |
596 |
|
|
rotation, bits, index, getIEEEBits); |
597 |
|
|
} |
598 |
|
|
|
599 |
|
|
/***************************************************************** |
600 |
|
|
* hilbert_box_vtx |
601 |
|
|
* |
602 |
|
|
* Determine the first or last vertex of a box to lie on a Hilbert curve |
603 |
|
|
* Inputs: |
604 |
|
|
* nDims: Number of coordinates. |
605 |
|
|
* nBytes: Number of bytes/coordinate. |
606 |
|
|
* nBits: Number of bits/coordinate. |
607 |
|
|
* findMin: Is it the least vertex sought? |
608 |
|
|
* coord1: Array of nDims nBytes-byte coordinates - one corner of box |
609 |
|
|
* coord2: Array of nDims nBytes-byte coordinates - opposite corner |
610 |
|
|
* Output: |
611 |
|
|
* c1 and c2 modified to refer to selected corner |
612 |
|
|
* value returned is log2 of size of largest power-of-two-aligned box that |
613 |
|
|
* contains the selected corner and no other corners |
614 |
|
|
* Assumptions: |
615 |
|
|
* nBits <= (sizeof bitmask_t) * (bits_per_byte) |
616 |
|
|
*/ |
617 |
|
|
|
618 |
|
|
|
619 |
|
|
static unsigned |
620 |
|
|
hilbert_box_vtx_work(unsigned nDims, unsigned nBytes, unsigned nBits, |
621 |
|
|
int findMin, |
622 |
|
|
unsigned max, unsigned y, |
623 |
|
|
char* c1, char* c2, |
624 |
|
|
unsigned rotation, |
625 |
|
|
bitmask_t bits, |
626 |
|
|
bitmask_t index, |
627 |
|
|
BitReader getBits) |
628 |
|
|
{ |
629 |
|
|
bitmask_t const one = 1; |
630 |
|
|
bitmask_t const ndOnes = ones(bitmask_t,nDims); |
631 |
|
|
bitmask_t const nd1Ones= ndOnes >> 1; |
632 |
|
|
bitmask_t bitsFolded = 0; |
633 |
|
|
|
634 |
|
|
while (y--) |
635 |
|
|
{ |
636 |
|
|
bitmask_t reflection = getBits(nDims, nBytes, c1, y); |
637 |
|
|
bitmask_t diff = reflection ^ getBits(nDims, nBytes, c2, y); |
638 |
|
|
if (diff) |
639 |
|
|
{ |
640 |
|
|
unsigned d; |
641 |
|
|
bitmask_t smear = rotateRight(diff, rotation, nDims) >> 1; |
642 |
|
|
bitmask_t digit = rotateRight(bits ^ reflection, rotation, nDims); |
643 |
|
|
for (d = 1; d < nDims; d *= 2) |
644 |
|
|
{ |
645 |
|
|
index ^= index >> d; |
646 |
|
|
digit ^= (digit >> d) &~ smear; |
647 |
|
|
smear |= smear >> d; |
648 |
|
|
} |
649 |
|
|
index &= 1; |
650 |
|
|
if ((index ^ y ^ findMin) & 1) |
651 |
|
|
digit ^= smear+1; |
652 |
|
|
digit = rotateLeft(digit, rotation, nDims) & diff; |
653 |
|
|
reflection ^= digit; |
654 |
|
|
|
655 |
|
|
for (d = 0; d < nDims; ++d) |
656 |
|
|
if (rdbit(diff, d)) |
657 |
|
|
{ |
658 |
|
|
int way = rdbit(digit, d); |
659 |
|
|
char* target = d*nBytes + (way? c1: c2); |
660 |
|
|
char* const source = 2*d*nBytes + c1 - target + c2; |
661 |
|
|
memcpy(target, source, nBytes); |
662 |
|
|
} |
663 |
|
|
|
664 |
|
|
bitsFolded |= diff; |
665 |
|
|
if (bitsFolded == ndOnes) |
666 |
|
|
return y; |
667 |
|
|
} |
668 |
|
|
|
669 |
|
|
bits ^= reflection; |
670 |
|
|
bits = rotateRight(bits, rotation, nDims); |
671 |
|
|
index ^= bits; |
672 |
|
|
reflection ^= one << rotation; |
673 |
|
|
adjust_rotation(rotation,nDims,bits); |
674 |
|
|
bits = reflection; |
675 |
|
|
} |
676 |
|
|
return y; |
677 |
|
|
} |
678 |
|
|
|
679 |
|
|
unsigned |
680 |
|
|
hilbert_box_vtx(unsigned nDims, unsigned nBytes, unsigned nBits, |
681 |
|
|
int findMin, void* c1, void* c2) |
682 |
|
|
{ |
683 |
|
|
bitmask_t const one = 1; |
684 |
|
|
bitmask_t bits = one << (nDims-1); |
685 |
|
|
return hilbert_box_vtx_work(nDims, nBytes, nBits, findMin, |
686 |
|
|
0, nBits, (char*)c1, (char*)c2, |
687 |
|
|
0, bits, bits, getIntBits); |
688 |
|
|
} |
689 |
|
|
|
690 |
|
|
unsigned |
691 |
|
|
hilbert_ieee_box_vtx(unsigned nDims, |
692 |
|
|
int findMin, double* c1, double* c2) |
693 |
|
|
{ |
694 |
|
|
unsigned rotation, max; |
695 |
|
|
bitmask_t bits, index; |
696 |
|
|
if (getIEEESignBits(nDims, c1) != getIEEESignBits(nDims, c2)) |
697 |
|
|
max = 2047; |
698 |
|
|
else |
699 |
|
|
{ |
700 |
|
|
unsigned max1 = getIEEEexptMax(nDims, c1); |
701 |
|
|
unsigned max2 = getIEEEexptMax(nDims, c2); |
702 |
|
|
max = (max1 > max2)? max1: max2; |
703 |
|
|
} |
704 |
|
|
|
705 |
|
|
getIEEEinitValues(c1, max+53, nDims, &rotation, &bits, &index); |
706 |
|
|
|
707 |
|
|
return hilbert_box_vtx_work(nDims, 8, 64, findMin, |
708 |
|
|
max, max+53, (char *)c1, (char *)c2, |
709 |
|
|
rotation, bits, index, getIEEEBits); |
710 |
|
|
} |
711 |
|
|
|
712 |
|
|
/***************************************************************** |
713 |
|
|
* hilbert_box_pt |
714 |
|
|
* |
715 |
|
|
* Determine the first or last point of a box to lie on a Hilbert curve |
716 |
|
|
* Inputs: |
717 |
|
|
* nDims: Number of coordinates. |
718 |
|
|
* nBytes: Number of bytes/coordinate. |
719 |
|
|
* nBits: Number of bits/coordinate. |
720 |
|
|
* findMin: Is it the least vertex sought? |
721 |
|
|
* coord1: Array of nDims nBytes-byte coordinates - one corner of box |
722 |
|
|
* coord2: Array of nDims nBytes-byte coordinates - opposite corner |
723 |
|
|
* Output: |
724 |
|
|
* c1 and c2 modified to refer to least point |
725 |
|
|
* Assumptions: |
726 |
|
|
* nBits <= (sizeof bitmask_t) * (bits_per_byte) |
727 |
|
|
*/ |
728 |
|
|
unsigned |
729 |
|
|
hilbert_box_pt_work(unsigned nDims, unsigned nBytes, unsigned nBits, |
730 |
|
|
int findMin, |
731 |
|
|
unsigned max, unsigned y, |
732 |
|
|
char* c1, char* c2, |
733 |
|
|
unsigned rotation, |
734 |
|
|
bitmask_t bits, |
735 |
|
|
bitmask_t index, |
736 |
|
|
BitReader getBits, |
737 |
|
|
BitWriter propogateBits) |
738 |
|
|
{ |
739 |
|
|
bitmask_t const one = 1; |
740 |
|
|
bitmask_t const nd1Ones = ones(bitmask_t,nDims) >> 1; |
741 |
|
|
bitmask_t fold1 = 0, fold2 = 0; |
742 |
|
|
unsigned smearSum = 0; |
743 |
|
|
|
744 |
|
|
while (y-- > max) |
745 |
|
|
{ |
746 |
|
|
bitmask_t reflection = getBits(nDims, nBytes, c1, y); |
747 |
|
|
bitmask_t diff = reflection ^ getBits(nDims, nBytes, c2, y); |
748 |
|
|
if (diff) |
749 |
|
|
{ |
750 |
|
|
bitmask_t smear = rotateRight(diff, rotation, nDims) >> 1; |
751 |
|
|
bitmask_t digit = rotateRight(bits ^ reflection, rotation, nDims); |
752 |
|
|
unsigned d; |
753 |
|
|
for (d = 1; d < nDims; d *= 2) |
754 |
|
|
{ |
755 |
|
|
index ^= index >> d; |
756 |
|
|
digit ^= (digit >> d) &~ smear; |
757 |
|
|
smear |= smear >> d; |
758 |
|
|
} |
759 |
|
|
smearSum += smear; |
760 |
|
|
index &= 1; |
761 |
|
|
if ((index ^ y ^ findMin) & 1) |
762 |
|
|
digit ^= smear+1; |
763 |
|
|
digit = rotateLeft(digit, rotation, nDims) & diff; |
764 |
|
|
reflection ^= digit; |
765 |
|
|
|
766 |
|
|
for (d = 0; d < nDims; ++d) |
767 |
|
|
if (rdbit(diff, d)) |
768 |
|
|
{ |
769 |
|
|
int way = rdbit(digit, d); |
770 |
|
|
char* c = way? c1: c2; |
771 |
|
|
bitmask_t fold = way? fold1: fold2; |
772 |
|
|
propogateBits(d, nBytes, c, y, rdbit(fold, d)); |
773 |
|
|
} |
774 |
|
|
diff ^= digit; |
775 |
|
|
fold1 |= digit; |
776 |
|
|
fold2 |= diff; |
777 |
|
|
} |
778 |
|
|
|
779 |
|
|
bits ^= reflection; |
780 |
|
|
bits = rotateRight(bits, rotation, nDims); |
781 |
|
|
index ^= bits; |
782 |
|
|
reflection ^= one << rotation; |
783 |
|
|
adjust_rotation(rotation,nDims,bits); |
784 |
|
|
bits = reflection; |
785 |
|
|
} |
786 |
|
|
return smearSum; |
787 |
|
|
} |
788 |
|
|
|
789 |
|
|
unsigned |
790 |
|
|
hilbert_box_pt(unsigned nDims, unsigned nBytes, unsigned nBits, |
791 |
|
|
int findMin, void* c1, void* c2) |
792 |
|
|
{ |
793 |
|
|
bitmask_t const one = 1; |
794 |
|
|
bitmask_t bits = one << (nDims-1); |
795 |
|
|
return hilbert_box_pt_work(nDims, nBytes, nBits, findMin, |
796 |
|
|
0, nBits, (char*)c1, (char*)c2, |
797 |
|
|
0, bits, bits, |
798 |
|
|
getIntBits, propogateIntBits); |
799 |
|
|
} |
800 |
|
|
|
801 |
|
|
unsigned |
802 |
|
|
hilbert_ieee_box_pt(unsigned nDims, |
803 |
|
|
int findMin, double* c1, double* c2) |
804 |
|
|
{ |
805 |
|
|
unsigned rotation, max; |
806 |
|
|
bitmask_t bits, index; |
807 |
|
|
bitmask_t c1Signs = getIEEESignBits(nDims, c1); |
808 |
|
|
bitmask_t c2Signs = getIEEESignBits(nDims, c2); |
809 |
|
|
if (c1Signs != c2Signs) |
810 |
|
|
{ |
811 |
|
|
rotation = 0; |
812 |
|
|
bits = (bitmask_t)1 << (nDims-1); |
813 |
|
|
index = 1; |
814 |
|
|
hilbert_box_pt_work(nDims, 8, 64, findMin, |
815 |
|
|
IEEErepBits-1, IEEErepBits, (char *)c1, (char *)c2, |
816 |
|
|
rotation, bits, index, |
817 |
|
|
getIEEEBits, propogateIEEEBits); |
818 |
|
|
} |
819 |
|
|
|
820 |
|
|
/* having put everything in the same orthant, start */ |
821 |
|
|
{ |
822 |
|
|
unsigned max1 = getIEEEexptMax(nDims, c1); |
823 |
|
|
unsigned max2 = getIEEEexptMax(nDims, c2); |
824 |
|
|
max = (max1 > max2)? max1: max2; |
825 |
|
|
} |
826 |
|
|
|
827 |
|
|
getIEEEinitValues(c1, max+53, nDims, &rotation, &bits, &index); |
828 |
|
|
|
829 |
|
|
return hilbert_box_pt_work(nDims, 8, 64, findMin, |
830 |
|
|
max, max+53, (char *)c1, (char *)c2, |
831 |
|
|
rotation, bits, index, |
832 |
|
|
getIEEEBits, propogateIEEEBits); |
833 |
|
|
} |
834 |
|
|
|
835 |
|
|
/***************************************************************** |
836 |
|
|
* hilbert_nextinbox |
837 |
|
|
* |
838 |
|
|
* Determine the first point of a box after or before a given point to lie on |
839 |
|
|
* a Hilbert curve |
840 |
|
|
* Inputs: |
841 |
|
|
* nDims: Number of coordinates. |
842 |
|
|
* nBytes: Number of bytes/coordinate. |
843 |
|
|
* nBits: Number of bits/coordinate. |
844 |
|
|
* findPrev: Is it a previous point that you want? |
845 |
|
|
* coord1: Array of nDims nBytes-byte coordinates - one corner of box |
846 |
|
|
* coord2: Array of nDims nBytes-byte coordinates - opposite corner |
847 |
|
|
* point: Array of nDims nBytes-byte coordinates - lower bound on point returned |
848 |
|
|
* |
849 |
|
|
* Output: |
850 |
|
|
if returns 1: |
851 |
|
|
* c1 and c2 modified to refer to least point after "point" in box |
852 |
|
|
else returns 0: |
853 |
|
|
arguments unchanged; "point" is beyond the last point of the box |
854 |
|
|
* Assumptions: |
855 |
|
|
* nBits <= (sizeof bitmask_t) * (bits_per_byte) |
856 |
|
|
*/ |
857 |
|
|
int |
858 |
|
|
hilbert_nextinbox(unsigned nDims, unsigned nBytes, unsigned nBits, |
859 |
|
|
int findPrev, void* c1V, void* c2V, void const* ptV) |
860 |
|
|
{ |
861 |
|
|
bitmask_t const one = 1; |
862 |
|
|
unsigned y = nBits; |
863 |
|
|
bitmask_t const ndOnes = ones(bitmask_t,nDims); |
864 |
|
|
bitmask_t const nd1Ones = ndOnes >> 1; |
865 |
|
|
unsigned rotation = 0; |
866 |
|
|
bitmask_t bits = 0; |
867 |
|
|
bitmask_t index = 0; |
868 |
|
|
bitmask_t fold1 = 0, fold2 = 0; |
869 |
|
|
bitmask_t valu1 = 0, valu2 = 0; |
870 |
|
|
unsigned p_y; |
871 |
|
|
bitmask_t p_separator = 0, p_firstSeparator; |
872 |
|
|
bitmask_t p_cornerdiff, p_reflection; |
873 |
|
|
bitmask_t p_fold1, p_fold2, p_valu1, p_valu2; |
874 |
|
|
|
875 |
|
|
char* c1 = (char*)c1V; |
876 |
|
|
char* c2 = (char*)c2V; |
877 |
|
|
char const* pt = (const char*)ptV; |
878 |
|
|
|
879 |
|
|
while (y-- > 0) |
880 |
|
|
{ |
881 |
|
|
bitmask_t reflection = getIntBits(nDims, nBytes, pt, y); |
882 |
|
|
bitmask_t diff = reflection ^ /* planes that separate box and point */ |
883 |
|
|
((getIntBits(nDims, nBytes, c1, y) &~ fold1) | valu1); |
884 |
|
|
|
885 |
|
|
if (diff) |
886 |
|
|
/* some coordinate planes separate point from box or |
887 |
|
|
dividing box or both; smear the bits of diff to reflect that |
888 |
|
|
after the first diff dimension, they might as well all be |
889 |
|
|
diffing; adjust the diff to reflect the fact that diffed |
890 |
|
|
dimensions don't matter. */ |
891 |
|
|
{ |
892 |
|
|
/* compute (the complement of) a "digit" in the integer index of this |
893 |
|
|
point */ |
894 |
|
|
bitmask_t cornerdiff = (diff ^ reflection) ^ /* separate box crnrs */ |
895 |
|
|
((getIntBits(nDims, nBytes, c2, y) &~ fold2) | valu2); |
896 |
|
|
bitmask_t separator = diff & ~cornerdiff; |
897 |
|
|
/* eventually, the most significant separating cutting plane */ |
898 |
|
|
bitmask_t firstSeparator; |
899 |
|
|
/* bits less significant than the msb of separator are irrelevant; |
900 |
|
|
for convenience, call them all separators too */ |
901 |
|
|
bitmask_t rotSep = rotateRight(separator, rotation, nDims); |
902 |
|
|
/* compute the (complement of the) digit of the hilbert code |
903 |
|
|
assoc with point */ |
904 |
|
|
bitmask_t digit = rotateRight(bits ^ reflection, rotation, nDims); |
905 |
|
|
unsigned d; |
906 |
|
|
for (d = 1; d < nDims; d *= 2) |
907 |
|
|
{ |
908 |
|
|
index ^= index >> d; |
909 |
|
|
digit ^= digit >> d; |
910 |
|
|
rotSep |= rotSep >> d; |
911 |
|
|
} |
912 |
|
|
index &= 1; |
913 |
|
|
digit &= rotSep; |
914 |
|
|
if ((index ^ y ^ findPrev) & 1) |
915 |
|
|
digit ^= rotSep; |
916 |
|
|
|
917 |
|
|
separator = rotateLeft(rotSep, rotation, nDims); |
918 |
|
|
rotSep -= rotSep >> 1; |
919 |
|
|
firstSeparator = rotateLeft(rotSep, rotation, nDims); |
920 |
|
|
/* forget about all the planes that split the box, except those that |
921 |
|
|
are more significant than the most significant separator. */ |
922 |
|
|
cornerdiff &= ~separator; |
923 |
|
|
|
924 |
|
|
if (cornerdiff && digit) |
925 |
|
|
/* some coordinate planes divide the box. Call the part of the |
926 |
|
|
box in the same orthant as the point "here" and the part of |
927 |
|
|
the box in the next (or previous) orthant "there". Remember |
928 |
|
|
what the "there" orthant of the box looks like in case it |
929 |
|
|
turns out that the curve doesn't reenter the box "here" after |
930 |
|
|
(before) passing thru point. Continue working with the |
931 |
|
|
"here" part. If there is no "there" there, skip it */ |
932 |
|
|
{ |
933 |
|
|
p_firstSeparator = digit & -digit; |
934 |
|
|
p_separator = 2*p_firstSeparator-1; |
935 |
|
|
p_separator = rotateLeft(p_separator, rotation, nDims); |
936 |
|
|
p_firstSeparator = rotateLeft(p_firstSeparator, rotation, nDims); |
937 |
|
|
p_cornerdiff = cornerdiff &~ (p_separator ^ p_firstSeparator); |
938 |
|
|
p_y = y; |
939 |
|
|
p_reflection = reflection ^ p_firstSeparator; |
940 |
|
|
p_fold1 = fold1; |
941 |
|
|
p_fold2 = fold2; |
942 |
|
|
p_valu1 = valu1; |
943 |
|
|
p_valu2 = valu2; |
944 |
|
|
} |
945 |
|
|
|
946 |
|
|
if (digit < rotSep) |
947 |
|
|
|
948 |
|
|
/* use next box */ |
949 |
|
|
{ |
950 |
|
|
if (!p_separator) return 0; /* no next point */ |
951 |
|
|
separator = p_separator; |
952 |
|
|
firstSeparator = p_firstSeparator; |
953 |
|
|
y = p_y; |
954 |
|
|
cornerdiff = p_cornerdiff; |
955 |
|
|
reflection = p_reflection; |
956 |
|
|
fold1 = p_fold1; |
957 |
|
|
fold2 = p_fold2; |
958 |
|
|
valu1 = p_valu1; |
959 |
|
|
valu2 = p_valu2; |
960 |
|
|
} |
961 |
|
|
|
962 |
|
|
if (cornerdiff) |
963 |
|
|
{ |
964 |
|
|
/* reduce currbox */ |
965 |
|
|
bitmask_t corner = diff & cornerdiff; |
966 |
|
|
cornerdiff ^= corner; |
967 |
|
|
fold1 |= corner; |
968 |
|
|
fold2 |= cornerdiff; |
969 |
|
|
valu1 |= ~reflection & corner; |
970 |
|
|
valu2 |= ~reflection & cornerdiff; |
971 |
|
|
} |
972 |
|
|
|
973 |
|
|
separator ^= firstSeparator; |
974 |
|
|
if (firstSeparator) |
975 |
|
|
/* we have completely separated the point from a part of the box |
976 |
|
|
ahead of it on the curve; almost done */ |
977 |
|
|
{ |
978 |
|
|
unsigned byteId = whichByte(nBytes,y); |
979 |
|
|
bitmask_t bthbit = one << y%8; |
980 |
|
|
for (d = 0; d < nDims; ++d) |
981 |
|
|
{ |
982 |
|
|
char lo1, lo2; |
983 |
|
|
char* cc1 = &c1[d*nBytes]; |
984 |
|
|
char* cc2 = &c2[d*nBytes]; |
985 |
|
|
char const* pnt = &pt[d*nBytes]; |
986 |
|
|
char hibits = -bthbit; |
987 |
|
|
char hipart = pnt[byteId] & hibits; |
988 |
|
|
memcpy(cc1, pnt, byteId); |
989 |
|
|
memcpy(cc2, pnt, byteId); |
990 |
|
|
|
991 |
|
|
if (rdbit(separator, d)) |
992 |
|
|
hibits ^= bthbit; |
993 |
|
|
if (rdbit(firstSeparator, d)) |
994 |
|
|
hipart ^= bthbit; |
995 |
|
|
|
996 |
|
|
if (rdbit(fold1, d)) |
997 |
|
|
{ |
998 |
|
|
lo1 = -rdbit(valu1, d); |
999 |
|
|
setBytes(cc1,byteId,nBytes,lo1); |
1000 |
|
|
} |
1001 |
|
|
else lo1 = cc1[byteId]; |
1002 |
|
|
cc1[byteId] = hipart | (lo1 &~ hibits); |
1003 |
|
|
|
1004 |
|
|
if (rdbit(fold2, d)) |
1005 |
|
|
{ |
1006 |
|
|
lo2 = -rdbit(valu2, d); |
1007 |
|
|
setBytes(cc2,byteId,nBytes,lo2); |
1008 |
|
|
} |
1009 |
|
|
else lo2 = cc2[byteId]; |
1010 |
|
|
cc2[byteId] = hipart | (lo2 &~ hibits); |
1011 |
|
|
} |
1012 |
|
|
|
1013 |
|
|
hilbert_box_pt(nDims, nBytes, nBits, !findPrev, c1V, c2V); |
1014 |
|
|
return 1; |
1015 |
|
|
} |
1016 |
|
|
} |
1017 |
|
|
|
1018 |
|
|
bits ^= reflection; |
1019 |
|
|
bits = rotateRight(bits, rotation, nDims); |
1020 |
|
|
index ^= bits; |
1021 |
|
|
reflection ^= one << rotation; |
1022 |
|
|
adjust_rotation(rotation,nDims,bits); |
1023 |
|
|
bits = reflection; |
1024 |
|
|
} |
1025 |
|
|
|
1026 |
|
|
/* point is in box */ |
1027 |
|
|
{ |
1028 |
|
|
unsigned d; |
1029 |
|
|
for (d = 0; d < nDims; ++d) |
1030 |
|
|
((char*)c1)[d] = ((char*)c2)[d] = ((char*)pt)[d]; |
1031 |
|
|
} |
1032 |
|
|
return 1; |
1033 |
|
|
} |
1034 |
|
|
|
1035 |
|
|
|
1036 |
|
|
|
1037 |
|
|
/***************************************************************** |
1038 |
|
|
* hilbert_incr |
1039 |
|
|
* |
1040 |
|
|
* Advance from one point to its successor on a Hilbert curve |
1041 |
|
|
* Inputs: |
1042 |
|
|
* nDims: Number of coordinates. |
1043 |
|
|
* nBits: Number of bits/coordinate. |
1044 |
|
|
* coord: Array of nDims nBits-bit coordinates. |
1045 |
|
|
* Output: |
1046 |
|
|
* coord: Next point on Hilbert curve |
1047 |
|
|
* Assumptions: |
1048 |
|
|
* nBits <= (sizeof bitmask_t) * (bits_per_byte) |
1049 |
|
|
*/ |
1050 |
|
|
|
1051 |
|
|
void |
1052 |
|
|
hilbert_incr(unsigned nDims, unsigned nBits, bitmask_t coord[]) |
1053 |
|
|
{ |
1054 |
|
|
bitmask_t const one = 1; |
1055 |
|
|
bitmask_t const ndOnes = ones(bitmask_t,nDims); |
1056 |
|
|
bitmask_t const nd1Ones= ndOnes >> 1; |
1057 |
|
|
unsigned b, d; |
1058 |
|
|
unsigned rotation = 0; |
1059 |
|
|
bitmask_t reflection = 0; |
1060 |
|
|
bitmask_t index = 0; |
1061 |
|
|
unsigned rb = nBits-1; |
1062 |
|
|
bitmask_t rd = ndOnes; |
1063 |
|
|
|
1064 |
|
|
for (b = nBits; b--;) |
1065 |
|
|
{ |
1066 |
|
|
bitmask_t bits = reflection; |
1067 |
|
|
reflection = 0; |
1068 |
|
|
for (d = 0; d < nDims; ++d) |
1069 |
|
|
reflection |= rdbit(coord[d], b) << d; |
1070 |
|
|
bits ^= reflection; |
1071 |
|
|
bits = rotateRight(bits, rotation, nDims); |
1072 |
|
|
index ^= bits; |
1073 |
|
|
for (d = 1; d < nDims; d *= 2) |
1074 |
|
|
index ^= index >> d; |
1075 |
|
|
if (index++ != ndOnes) |
1076 |
|
|
{ |
1077 |
|
|
rb = b; |
1078 |
|
|
rd = index & -index; |
1079 |
|
|
rd = rotateLeft(rd, rotation, nDims); |
1080 |
|
|
|
1081 |
|
|
} |
1082 |
|
|
index &= 1; |
1083 |
|
|
index <<= nDims-1; |
1084 |
|
|
|
1085 |
|
|
reflection ^= one << rotation; |
1086 |
|
|
adjust_rotation(rotation,nDims,bits); |
1087 |
|
|
} |
1088 |
|
|
for (d = 0; !rdbit(rd, d); ++d) {} |
1089 |
|
|
coord[d] ^= (2 << rb) - 1; |
1090 |
|
|
} |
1091 |
|
|
|
1092 |
|
|
|
1093 |
|
|
/* LICENSE |
1094 |
|
|
* |
1095 |
|
|
* This software is copyrighted by Rice University. It may be freely copied, |
1096 |
|
|
* modified, and redistributed, provided that the copyright notice is |
1097 |
|
|
* preserved on all copies. |
1098 |
|
|
* |
1099 |
|
|
* There is no warranty or other guarantee of fitness for this software, |
1100 |
|
|
* it is provided solely "as is". Bug reports or fixes may be sent |
1101 |
|
|
* to the author, who may or may not act on them as he desires. |
1102 |
|
|
* |
1103 |
|
|
* You may include this software in a program or other software product, |
1104 |
|
|
* but must display the notice: |
1105 |
|
|
* |
1106 |
|
|
* Hilbert Curve implementation copyright 1998, Rice University |
1107 |
|
|
* |
1108 |
|
|
* in any place where the end-user would see your own copyright. |
1109 |
|
|
* |
1110 |
|
|
* If you modify this software, you should include a notice giving the |
1111 |
|
|
* name of the person performing the modification, the date of modification, |
1112 |
|
|
* and the reason for such modification. |
1113 |
|
|
*/ |
1114 |
|
|
|
1115 |
|
|
|
1116 |
|
|
|
1117 |
|
|
/* Revision history: |
1118 |
|
|
|
1119 |
|
|
July 1998: Initial release |
1120 |
|
|
|
1121 |
|
|
Sept 1998: Second release |
1122 |
|
|
|
1123 |
|
|
Dec 1998: Fixed bug in hilbert_c2i that allowed a shift by number of bits in |
1124 |
|
|
bitmask to vaporize index, in last bit of the function. Implemented |
1125 |
|
|
hilbert_incr. |
1126 |
|
|
|
1127 |
|
|
August 1999: Added argument to hilbert_nextinbox so that you can, optionally, |
1128 |
|
|
find the previous point along the curve to intersect the box, rather than the |
1129 |
|
|
next point. |
1130 |
|
|
|
1131 |
|
|
Nov 1999: Defined fast bit-transpose function (fast, at least, if the number |
1132 |
|
|
of bits is large), and reimplemented i2c and c2i in terms of it. Collapsed |
1133 |
|
|
loops in hilbert_cmp, with the intention of reusing the cmp code to compare |
1134 |
|
|
more general bitstreams. |
1135 |
|
|
|
1136 |
|
|
Feb 2000: Implemented almost all the floating point versions of cmp, etc, so |
1137 |
|
|
that coordinates expressed in terms of double-precision IEEE floating point |
1138 |
|
|
can be ordered. Still have to do next-in-box, though. |
1139 |
|
|
|
1140 |
|
|
Oct 2001: Learned that some arbitrary coding choices caused some routines |
1141 |
|
|
to fail in one dimension, and changed those choices. |
1142 |
|
|
|
1143 |
|
|
version 2001-10-20-05:34 |
1144 |
|
|
|
1145 |
|
|
*/ |
1146 |
|
|
|
1147 |
|
|
/* What remains is test code that won't be compiled unless you define the |
1148 |
|
|
TEST_HILBERT preprocessor symbol */ |
1149 |
|
|
|
1150 |
|
|
#ifdef TEST_HILBERT |
1151 |
|
|
#include <stdio.h> |
1152 |
|
|
#define abs(x) (((x)>=0)?(x):(-(x))) |
1153 |
|
|
|
1154 |
|
|
int main() |
1155 |
|
|
{ |
1156 |
|
|
#define maxDim (8*sizeof(bitmask_t)) |
1157 |
|
|
bitmask_t coord[maxDim], coordPrev[maxDim]; |
1158 |
|
|
unsigned nDims, nBits, nPrints, orderCheck, i; |
1159 |
|
|
bitmask_t r, r1; |
1160 |
|
|
|
1161 |
|
|
for (;;) |
1162 |
|
|
{ |
1163 |
|
|
printf( "Enter nDims, nBits, nPrints, orderCheck: " ); |
1164 |
|
|
scanf( "%d", &nDims); |
1165 |
|
|
if ( nDims == 0 ) |
1166 |
|
|
break; |
1167 |
|
|
scanf( "%d%d%d", &nBits, &nPrints, &orderCheck); |
1168 |
|
|
while ( (i = getchar()) != '\n' && i != EOF ) |
1169 |
|
|
; |
1170 |
|
|
if ( i == EOF ) |
1171 |
|
|
break; |
1172 |
|
|
|
1173 |
|
|
if (nDims*nBits > 8*sizeof(r)) |
1174 |
|
|
{ |
1175 |
|
|
printf("Product of nDims and nBits not exceed %d.\n", 8*sizeof(r)); |
1176 |
|
|
break; |
1177 |
|
|
} |
1178 |
|
|
|
1179 |
|
|
if (nBits == 0) |
1180 |
|
|
{ |
1181 |
|
|
printf("nBits must be positive.\n"); |
1182 |
|
|
break; |
1183 |
|
|
} |
1184 |
|
|
|
1185 |
|
|
if (nPrints > (1ULL << (nDims*nBits))) |
1186 |
|
|
nPrints = 1ULL << (nDims*nBits); |
1187 |
|
|
|
1188 |
|
|
for (r = 0; r < nPrints; ++r) |
1189 |
|
|
{ |
1190 |
|
|
bitmask_t coord1[maxDim]; |
1191 |
|
|
int miscount = 0; |
1192 |
|
|
hilbert_i2c( nDims, nBits, r, coord ); |
1193 |
|
|
printf("%d: ", (unsigned)r); |
1194 |
|
|
for (i = 0; i < nDims; ++i) |
1195 |
|
|
{ |
1196 |
|
|
int diff = (int)(coord[i] - coordPrev[i]); |
1197 |
|
|
miscount += abs(diff); |
1198 |
|
|
coordPrev[i] = coord[i]; |
1199 |
|
|
printf(" %d", (unsigned)coord[i]); |
1200 |
|
|
} |
1201 |
|
|
if (r > 0 && miscount != 1) |
1202 |
|
|
printf(".....error"); |
1203 |
|
|
printf("\n"); |
1204 |
|
|
r1 = hilbert_c2i( nDims, nBits, coord ); |
1205 |
|
|
if ( r != r1 ) |
1206 |
|
|
printf( "r = 0x%x; r1 = 0x%x\n", (unsigned)r, (unsigned)r1); |
1207 |
|
|
for (i = 0; i < nDims; ++i) |
1208 |
|
|
coord[i] = coordPrev[i]; |
1209 |
|
|
|
1210 |
|
|
if (! orderCheck) |
1211 |
|
|
continue; |
1212 |
|
|
|
1213 |
|
|
for (r1 = 0; r1 < r; ++r1 ) |
1214 |
|
|
{ |
1215 |
|
|
unsigned ans; |
1216 |
|
|
hilbert_i2c( nDims, nBits, r1, coord1 ); |
1217 |
|
|
ans = hilbert_cmp( nDims, sizeof(coord[0]), nBits, coord, coord1); |
1218 |
|
|
if (ans != 1) |
1219 |
|
|
{ |
1220 |
|
|
int width = (nDims*nBits + 3) / 4; |
1221 |
|
|
printf( "cmp r = 0x%0*x; r1 = 0x%0*x, ans = %2d\n", |
1222 |
|
|
width, (unsigned)r, |
1223 |
|
|
width, (unsigned)r1, ans ); |
1224 |
|
|
} |
1225 |
|
|
} |
1226 |
|
|
hilbert_i2c( nDims, nBits, r1, coord1 ); |
1227 |
|
|
if (hilbert_cmp( nDims, sizeof(coord[0]), nBits, coord, coord1) != 0) |
1228 |
|
|
printf( "cmp r = 0x%0*x; r1 = 0x%0*x\n", (nDims*nBits+3)/4, (unsigned)r, |
1229 |
|
|
(nDims*nBits+3)/4, (unsigned)r1 ); |
1230 |
|
|
|
1231 |
|
|
} |
1232 |
|
|
} |
1233 |
|
|
return 0; |
1234 |
|
|
} |
1235 |
|
|
|
1236 |
|
|
#endif |
1237 |
|
|
|
1238 |
|
|
#ifdef TEST_IEEE |
1239 |
|
|
#include <stdio.h> |
1240 |
|
|
#include <stdlib.h> |
1241 |
|
|
#include <math.h> |
1242 |
|
|
|
1243 |
|
|
int cmp(const void* xv, const void* yv) |
1244 |
|
|
{ |
1245 |
|
|
double const* x = xv; |
1246 |
|
|
double const* y = yv; |
1247 |
|
|
/* return hilbert_cmp(2, 8, 64, x, y); */ |
1248 |
|
|
return hilbert_ieee_cmp(2, x, y); |
1249 |
|
|
} |
1250 |
|
|
|
1251 |
|
|
int main() |
1252 |
|
|
{ |
1253 |
|
|
double *a; |
1254 |
|
|
unsigned i; |
1255 |
|
|
unsigned n; |
1256 |
|
|
printf("How many points? "); |
1257 |
|
|
scanf("%d", &n); |
1258 |
|
|
a = (double*) malloc(2*n*sizeof(double)); |
1259 |
|
|
for (i = 0; i < n; ++i) |
1260 |
|
|
a[2*i] = drand48()-0.5, a[2*i+1] = drand48()-0.5; |
1261 |
|
|
|
1262 |
|
|
qsort(a, n, 2*sizeof(double), cmp); |
1263 |
|
|
|
1264 |
|
|
for (i = 0; i < n; ++i) |
1265 |
|
|
printf("%8g %8g\n", a[2*i], a[2*i+1]); |
1266 |
|
|
free(a); |
1267 |
|
|
return 0; |
1268 |
|
|
} |
1269 |
|
|
|
1270 |
|
|
#endif |
1271 |
|
|
|
1272 |
|
|
#ifdef TEST_CMP |
1273 |
|
|
#include <stdio.h> |
1274 |
|
|
|
1275 |
|
|
#define maxDim (8*sizeof(bitmask_t)) |
1276 |
|
|
int main() |
1277 |
|
|
{ |
1278 |
|
|
double coord[maxDim]; |
1279 |
|
|
unsigned nDims, i, k; |
1280 |
|
|
|
1281 |
|
|
printf( "Enter nDims: " ); |
1282 |
|
|
scanf( "%d", &nDims); |
1283 |
|
|
if ( nDims == 0 ) |
1284 |
|
|
return 0; |
1285 |
|
|
while ( (i = getchar()) != '\n' && i != EOF ) |
1286 |
|
|
; |
1287 |
|
|
if ( i == EOF ) |
1288 |
|
|
return 0; |
1289 |
|
|
|
1290 |
|
|
for (k = 0; k < (1<<nDims); ++k) |
1291 |
|
|
{ |
1292 |
|
|
printf("Orth %2d\n", k); |
1293 |
|
|
for (i = 0; i < nDims; ++i) |
1294 |
|
|
coord[i] = ((k>>i)&1)? -1.: 1.; |
1295 |
|
|
|
1296 |
|
|
|
1297 |
|
|
hilbert_ieee_cmp( nDims, coord, coord); |
1298 |
|
|
} |
1299 |
|
|
return 0; |
1300 |
|
|
} |
1301 |
|
|
|
1302 |
|
|
#endif |
1303 |
|
|
|
1304 |
|
|
#ifdef TEST_VTX |
1305 |
|
|
#include <stdio.h> |
1306 |
|
|
#include <stdlib.h> |
1307 |
|
|
|
1308 |
|
|
#define maxDim (8*sizeof(bitmask_t)) |
1309 |
|
|
|
1310 |
|
|
unsigned g_nDims; |
1311 |
|
|
|
1312 |
|
|
int cmp(void const* c1p, void const* c2p) |
1313 |
|
|
{ |
1314 |
|
|
return hilbert_cmp(g_nDims, sizeof(unsigned), 8*sizeof(unsigned), c1p, c2p); |
1315 |
|
|
} |
1316 |
|
|
|
1317 |
|
|
int main() |
1318 |
|
|
{ |
1319 |
|
|
unsigned corner0[maxDim], corner1[maxDim]; |
1320 |
|
|
unsigned cornerlo[maxDim], cornerhi[maxDim], work[maxDim]; |
1321 |
|
|
typedef unsigned array_t[maxDim]; |
1322 |
|
|
array_t* array; |
1323 |
|
|
|
1324 |
|
|
unsigned nDims, i, k; |
1325 |
|
|
|
1326 |
|
|
printf( "Enter nDims: " ); |
1327 |
|
|
scanf( "%d", &nDims); |
1328 |
|
|
if ( nDims == 0 ) |
1329 |
|
|
return 0; |
1330 |
|
|
while ( (i = getchar()) != '\n' && i != EOF ) |
1331 |
|
|
; |
1332 |
|
|
if ( i == EOF ) |
1333 |
|
|
return 0; |
1334 |
|
|
|
1335 |
|
|
printf("Enter one corner (%d coordinates): ", nDims); |
1336 |
|
|
for (k = 0; k < nDims; ++k) |
1337 |
|
|
scanf("%d", &corner0[k]); |
1338 |
|
|
|
1339 |
|
|
printf("Enter other corner (%d coordinates): ", nDims); |
1340 |
|
|
for (k = 0; k < nDims; ++k) |
1341 |
|
|
scanf("%d", &corner1[k]); |
1342 |
|
|
|
1343 |
|
|
|
1344 |
|
|
/* find first corner */ |
1345 |
|
|
for (k = 0; k < nDims; ++k) |
1346 |
|
|
{ |
1347 |
|
|
cornerlo[k] = corner0[k]; |
1348 |
|
|
work[k] = corner1[k]; |
1349 |
|
|
} |
1350 |
|
|
|
1351 |
|
|
hilbert_box_vtx(nDims, sizeof(unsigned), 8*sizeof(unsigned), |
1352 |
|
|
1, cornerlo, work); |
1353 |
|
|
printf("Predicted lo corner: "); |
1354 |
|
|
for (k = 0; k < nDims; ++k) |
1355 |
|
|
printf("%4u", cornerlo[k]); |
1356 |
|
|
printf("\n"); |
1357 |
|
|
|
1358 |
|
|
|
1359 |
|
|
/* find last corner */ |
1360 |
|
|
for (k = 0; k < nDims; ++k) |
1361 |
|
|
{ |
1362 |
|
|
work[k] = corner0[k]; |
1363 |
|
|
cornerhi[k] = corner1[k]; |
1364 |
|
|
} |
1365 |
|
|
|
1366 |
|
|
hilbert_box_vtx(nDims, sizeof(unsigned), 8*sizeof(unsigned), |
1367 |
|
|
0, work, cornerhi); |
1368 |
|
|
printf("Predicted hi corner: "); |
1369 |
|
|
for (k = 0; k < nDims; ++k) |
1370 |
|
|
printf("%4u", cornerhi[k]); |
1371 |
|
|
printf("\n"); |
1372 |
|
|
|
1373 |
|
|
array = (array_t*) malloc(maxDim*sizeof(unsigned) << nDims); |
1374 |
|
|
for (k = 0; k < (1<<nDims); ++k) |
1375 |
|
|
{ |
1376 |
|
|
unsigned j; |
1377 |
|
|
unsigned* eltk = &array[k][0]; |
1378 |
|
|
for (j = 0; j < nDims; ++j) |
1379 |
|
|
{ |
1380 |
|
|
unsigned* src = ((k>>j)&1)? corner1: corner0; |
1381 |
|
|
eltk[j] = src[j]; |
1382 |
|
|
} |
1383 |
|
|
} |
1384 |
|
|
|
1385 |
|
|
g_nDims = nDims; |
1386 |
|
|
qsort(array, (1<<nDims), maxDim*sizeof(unsigned), cmp); |
1387 |
|
|
|
1388 |
|
|
printf("Result of sort\n"); |
1389 |
|
|
for (k = 0; k < (1<<nDims); k += (1 << nDims) - 1) |
1390 |
|
|
{ |
1391 |
|
|
unsigned j; |
1392 |
|
|
unsigned* eltk = &array[k][0]; |
1393 |
|
|
for (j = 0; j < nDims; ++j) |
1394 |
|
|
printf("%4u", eltk[j]); |
1395 |
|
|
printf("\n"); |
1396 |
|
|
} |
1397 |
|
|
free((char*)array); |
1398 |
|
|
return 0; |
1399 |
|
|
} |
1400 |
|
|
|
1401 |
|
|
#endif |
1402 |
|
|
|
1403 |
|
|
#ifdef TEST_IEEE_VTX |
1404 |
|
|
#include <stdio.h> |
1405 |
|
|
#include <stdlib.h> |
1406 |
|
|
#include <assert.h> |
1407 |
|
|
|
1408 |
|
|
#define maxDim (8*sizeof(bitmask_t)) |
1409 |
|
|
typedef double key_t; |
1410 |
|
|
|
1411 |
|
|
unsigned g_nDims; |
1412 |
|
|
|
1413 |
|
|
int cmp(void const* c1p, void const* c2p) |
1414 |
|
|
{ |
1415 |
|
|
return hilbert_ieee_cmp(g_nDims, c1p, c2p); |
1416 |
|
|
} |
1417 |
|
|
|
1418 |
|
|
int main() |
1419 |
|
|
{ |
1420 |
|
|
key_t corner0[maxDim], corner1[maxDim]; |
1421 |
|
|
key_t cornerlo[maxDim], cornerhi[maxDim], work[maxDim]; |
1422 |
|
|
typedef key_t array_t[maxDim]; |
1423 |
|
|
array_t* array; |
1424 |
|
|
|
1425 |
|
|
unsigned nDims, i, k; |
1426 |
|
|
|
1427 |
|
|
printf( "Enter nDims: " ); |
1428 |
|
|
scanf( "%d", &nDims); |
1429 |
|
|
if ( nDims == 0 ) |
1430 |
|
|
return 0; |
1431 |
|
|
|
1432 |
|
|
for (i = 0; i < 10000; ++i) |
1433 |
|
|
{ |
1434 |
|
|
for (k = 0; k < nDims; ++k) |
1435 |
|
|
{ |
1436 |
|
|
corner0[k] = 2.*drand48() - 1.; |
1437 |
|
|
corner1[k] = 2.*drand48() - 1.; |
1438 |
|
|
} |
1439 |
|
|
|
1440 |
|
|
/* find first corner */ |
1441 |
|
|
for (k = 0; k < nDims; ++k) |
1442 |
|
|
{ |
1443 |
|
|
cornerlo[k] = corner0[k]; |
1444 |
|
|
work[k] = corner1[k]; |
1445 |
|
|
} |
1446 |
|
|
|
1447 |
|
|
hilbert_ieee_box_vtx(nDims, 1, cornerlo, work); |
1448 |
|
|
|
1449 |
|
|
/* find last corner */ |
1450 |
|
|
for (k = 0; k < nDims; ++k) |
1451 |
|
|
{ |
1452 |
|
|
work[k] = corner0[k]; |
1453 |
|
|
cornerhi[k] = corner1[k]; |
1454 |
|
|
} |
1455 |
|
|
|
1456 |
|
|
hilbert_ieee_box_vtx(nDims, 0, work, cornerhi); |
1457 |
|
|
|
1458 |
|
|
array = (array_t*) malloc(maxDim*sizeof(key_t) << nDims); |
1459 |
|
|
for (k = 0; k < (1<<nDims); ++k) |
1460 |
|
|
{ |
1461 |
|
|
unsigned j; |
1462 |
|
|
key_t* eltk = &array[k][0]; |
1463 |
|
|
for (j = 0; j < nDims; ++j) |
1464 |
|
|
{ |
1465 |
|
|
key_t* src = ((k>>j)&1)? corner1: corner0; |
1466 |
|
|
eltk[j] = src[j]; |
1467 |
|
|
} |
1468 |
|
|
} |
1469 |
|
|
|
1470 |
|
|
g_nDims = nDims; |
1471 |
|
|
qsort(array, (1<<nDims), maxDim*sizeof(key_t), cmp); |
1472 |
|
|
|
1473 |
|
|
for (k = 0; k < (1<<nDims); k += (1 << nDims) - 1) |
1474 |
|
|
{ |
1475 |
|
|
unsigned j; |
1476 |
|
|
int mismatch = 0; |
1477 |
|
|
key_t* eltk = &array[k][0]; |
1478 |
|
|
for (j = 0; j < nDims & !mismatch; ++j) |
1479 |
|
|
{ |
1480 |
|
|
mismatch = (eltk[j] != ((k==0)? cornerlo: cornerhi)[j]); |
1481 |
|
|
} |
1482 |
|
|
assert (!mismatch); |
1483 |
|
|
} |
1484 |
|
|
free((char*)array); |
1485 |
|
|
} |
1486 |
|
|
return 0; |
1487 |
|
|
} |
1488 |
|
|
|
1489 |
|
|
#endif |
1490 |
|
|
|
1491 |
|
|
#ifdef TEST_PT |
1492 |
|
|
#include <stdio.h> |
1493 |
|
|
#include <stdlib.h> |
1494 |
|
|
|
1495 |
|
|
#define maxDim (8*sizeof(bitmask_t)) |
1496 |
|
|
|
1497 |
|
|
unsigned g_nDims; |
1498 |
|
|
|
1499 |
|
|
int cmp(void const* c1p, void const* c2p) |
1500 |
|
|
{ |
1501 |
|
|
return hilbert_cmp(g_nDims, sizeof(unsigned), 8*sizeof(unsigned), c1p, c2p); |
1502 |
|
|
} |
1503 |
|
|
|
1504 |
|
|
int main() |
1505 |
|
|
{ |
1506 |
|
|
unsigned point0[maxDim], point1[maxDim]; |
1507 |
|
|
unsigned pointlo[maxDim], pointhi[maxDim], work[maxDim]; |
1508 |
|
|
typedef unsigned array_t[maxDim]; |
1509 |
|
|
array_t* array; |
1510 |
|
|
|
1511 |
|
|
unsigned nDims, i, k, outvolume = 1, involume = 1; |
1512 |
|
|
unsigned nextItem; |
1513 |
|
|
|
1514 |
|
|
printf( "Enter nDims: " ); |
1515 |
|
|
scanf( "%d", &nDims); |
1516 |
|
|
if ( nDims == 0 ) |
1517 |
|
|
return 0; |
1518 |
|
|
while ( (i = getchar()) != '\n' && i != EOF ) |
1519 |
|
|
; |
1520 |
|
|
if ( i == EOF ) |
1521 |
|
|
return 0; |
1522 |
|
|
|
1523 |
|
|
printf("Enter one point (%d coordinates): ", nDims); |
1524 |
|
|
for (k = 0; k < nDims; ++k) |
1525 |
|
|
scanf("%d", &point0[k]); |
1526 |
|
|
|
1527 |
|
|
printf("Enter other point (%d coordinates, strictly greater): ", nDims); |
1528 |
|
|
for (k = 0; k < nDims; ++k) |
1529 |
|
|
{ |
1530 |
|
|
unsigned diff; |
1531 |
|
|
scanf("%d", &point1[k]); |
1532 |
|
|
diff = point1[k] - point0[k]; |
1533 |
|
|
outvolume *= diff + 1; |
1534 |
|
|
involume *= diff - 1; |
1535 |
|
|
} |
1536 |
|
|
|
1537 |
|
|
|
1538 |
|
|
/* find first point */ |
1539 |
|
|
for (k = 0; k < nDims; ++k) |
1540 |
|
|
{ |
1541 |
|
|
pointlo[k] = point0[k]; |
1542 |
|
|
work[k] = point1[k]; |
1543 |
|
|
} |
1544 |
|
|
|
1545 |
|
|
hilbert_box_pt(nDims, sizeof(unsigned), 8*sizeof(unsigned), |
1546 |
|
|
1, pointlo, work); |
1547 |
|
|
printf("Predicted lo point: "); |
1548 |
|
|
for (k = 0; k < nDims; ++k) |
1549 |
|
|
printf("%4u", pointlo[k]); |
1550 |
|
|
printf("\n"); |
1551 |
|
|
|
1552 |
|
|
|
1553 |
|
|
/* find last point */ |
1554 |
|
|
for (k = 0; k < nDims; ++k) |
1555 |
|
|
{ |
1556 |
|
|
work[k] = point0[k]; |
1557 |
|
|
pointhi[k] = point1[k]; |
1558 |
|
|
} |
1559 |
|
|
|
1560 |
|
|
hilbert_box_pt(nDims, sizeof(unsigned), 8*sizeof(unsigned), |
1561 |
|
|
0, work, pointhi); |
1562 |
|
|
printf("Predicted hi point: "); |
1563 |
|
|
for (k = 0; k < nDims; ++k) |
1564 |
|
|
printf("%4u", pointhi[k]); |
1565 |
|
|
printf("\n"); |
1566 |
|
|
|
1567 |
|
|
|
1568 |
|
|
|
1569 |
|
|
array = (array_t*) malloc(maxDim*sizeof(unsigned) * (outvolume-involume)); |
1570 |
|
|
if (array == 0) |
1571 |
|
|
{ |
1572 |
|
|
fprintf(stderr, "Out of memory.\n"); |
1573 |
|
|
exit(-1); |
1574 |
|
|
} |
1575 |
|
|
nextItem = 0; |
1576 |
|
|
for (k = 0; k < outvolume; ++k) |
1577 |
|
|
{ |
1578 |
|
|
unsigned kk = k; |
1579 |
|
|
unsigned j; |
1580 |
|
|
unsigned* eltk = &array[nextItem][0]; |
1581 |
|
|
int boundary = 0; |
1582 |
|
|
|
1583 |
|
|
for (j = 0; j < nDims; ++j) |
1584 |
|
|
{ |
1585 |
|
|
unsigned diff1 = point1[j] - point0[j] + 1; |
1586 |
|
|
unsigned pos = point0[j] + (kk % diff1); |
1587 |
|
|
boundary |= (point0[j] == pos || pos == point1[j]); |
1588 |
|
|
eltk[j] = pos; |
1589 |
|
|
kk /= diff1; |
1590 |
|
|
} |
1591 |
|
|
if (boundary) |
1592 |
|
|
++nextItem; |
1593 |
|
|
} |
1594 |
|
|
|
1595 |
|
|
g_nDims = nDims; |
1596 |
|
|
qsort(array, outvolume-involume, maxDim*sizeof(unsigned), cmp); |
1597 |
|
|
|
1598 |
|
|
printf("Result of sort\n"); |
1599 |
|
|
for (k = 0; k < outvolume-involume; k += outvolume-involume-1) |
1600 |
|
|
{ |
1601 |
|
|
unsigned j; |
1602 |
|
|
unsigned* eltk = &array[k][0]; |
1603 |
|
|
for (j = 0; j < nDims; ++j) |
1604 |
|
|
printf("%4u", eltk[j]); |
1605 |
|
|
printf("\n"); |
1606 |
|
|
} |
1607 |
|
|
free((char*)array); |
1608 |
|
|
return 0; |
1609 |
|
|
} |
1610 |
|
|
|
1611 |
|
|
#endif |
1612 |
|
|
|
1613 |
|
|
#ifdef TEST_IEEE_PT |
1614 |
|
|
#include <stdio.h> |
1615 |
|
|
#include <stdlib.h> |
1616 |
|
|
|
1617 |
|
|
#define maxDim (8*sizeof(bitmask_t)) |
1618 |
|
|
|
1619 |
|
|
int main() |
1620 |
|
|
{ |
1621 |
|
|
double point0[maxDim], point1[maxDim]; |
1622 |
|
|
double pointlo[maxDim], pointhi[maxDim], work[maxDim]; |
1623 |
|
|
|
1624 |
|
|
unsigned nDims, k, i; |
1625 |
|
|
|
1626 |
|
|
printf( "Enter nDims: " ); |
1627 |
|
|
scanf( "%d", &nDims); |
1628 |
|
|
if ( nDims == 0 ) |
1629 |
|
|
return 0; |
1630 |
|
|
while ( (i = getchar()) != '\n' && i != EOF ) |
1631 |
|
|
; |
1632 |
|
|
if ( i == EOF ) |
1633 |
|
|
return 0; |
1634 |
|
|
|
1635 |
|
|
printf("Enter one point (%d coordinates): ", nDims); |
1636 |
|
|
for (k = 0; k < nDims; ++k) |
1637 |
|
|
scanf("%lf", &point0[k]); |
1638 |
|
|
|
1639 |
|
|
printf("Enter other point (%d coordinates, strictly greater): ", nDims); |
1640 |
|
|
for (k = 0; k < nDims; ++k) |
1641 |
|
|
scanf("%lf", &point1[k]); |
1642 |
|
|
|
1643 |
|
|
/* find last point */ |
1644 |
|
|
for (k = 0; k < nDims; ++k) |
1645 |
|
|
{ |
1646 |
|
|
work[k] = point0[k]; |
1647 |
|
|
pointhi[k] = point1[k]; |
1648 |
|
|
} |
1649 |
|
|
|
1650 |
|
|
hilbert_ieee_box_pt(nDims, 0, work, pointhi); |
1651 |
|
|
printf("Predicted hi point: "); |
1652 |
|
|
for (k = 0; k < nDims; ++k) |
1653 |
|
|
printf("%10lg", pointhi[k]); |
1654 |
|
|
printf("\n"); |
1655 |
|
|
|
1656 |
|
|
/* find first point */ |
1657 |
|
|
for (k = 0; k < nDims; ++k) |
1658 |
|
|
{ |
1659 |
|
|
pointlo[k] = point0[k]; |
1660 |
|
|
work[k] = point1[k]; |
1661 |
|
|
} |
1662 |
|
|
|
1663 |
|
|
hilbert_ieee_box_pt(nDims, 1, pointlo, work); |
1664 |
|
|
printf("Predicted lo point: "); |
1665 |
|
|
for (k = 0; k < nDims; ++k) |
1666 |
|
|
printf("%10lg", pointlo[k]); |
1667 |
|
|
printf("\n"); |
1668 |
|
|
|
1669 |
|
|
/* validate by sorting random boundary points */ |
1670 |
|
|
#define nPts 1000000 |
1671 |
|
|
assert(hilbert_ieee_cmp(nDims, pointlo, pointhi) < 0); |
1672 |
|
|
for (i = 0; i < nPts; ++i) |
1673 |
|
|
{ |
1674 |
|
|
double pt1[maxDim], pt2[maxDim]; |
1675 |
|
|
for (k = 0; k < nDims; ++k) |
1676 |
|
|
{ |
1677 |
|
|
if (i % nDims == k) |
1678 |
|
|
pt1[k] = point0[k]; |
1679 |
|
|
else |
1680 |
|
|
pt1[k] = point0[k] + drand48()*(point1[k]-point0[k]); |
1681 |
|
|
} |
1682 |
|
|
for (k = 0; k < nDims; ++k) |
1683 |
|
|
{ |
1684 |
|
|
if (i % nDims == k) |
1685 |
|
|
pt2[k] = point1[k]; |
1686 |
|
|
else |
1687 |
|
|
pt2[k] = point0[k] + drand48()*(point1[k]-point0[k]); |
1688 |
|
|
} |
1689 |
|
|
if (hilbert_ieee_cmp(nDims, pt1, pt2) < 0) |
1690 |
|
|
{ |
1691 |
|
|
if (hilbert_ieee_cmp(nDims, pt1, pointlo) < 0) |
1692 |
|
|
memcpy(pointlo, pt1, maxDim*sizeof(double)); |
1693 |
|
|
if (hilbert_ieee_cmp(nDims, pointhi, pt2) < 0) |
1694 |
|
|
memcpy(pointhi, pt2, maxDim*sizeof(double)); |
1695 |
|
|
} |
1696 |
|
|
else |
1697 |
|
|
{ |
1698 |
|
|
if (hilbert_ieee_cmp(nDims, pt2, pointlo) < 0) |
1699 |
|
|
memcpy(pointlo, pt2, maxDim*sizeof(double)); |
1700 |
|
|
if (hilbert_ieee_cmp(nDims, pointhi, pt1) < 0) |
1701 |
|
|
memcpy(pointhi, pt1, maxDim*sizeof(double)); |
1702 |
|
|
} |
1703 |
|
|
} |
1704 |
|
|
|
1705 |
|
|
printf("Sorted hi and lo:\n"); |
1706 |
|
|
for (k = 0; k < nDims; ++k) |
1707 |
|
|
printf("%10lg", pointhi[k]); |
1708 |
|
|
printf("\n"); |
1709 |
|
|
for (k = 0; k < nDims; ++k) |
1710 |
|
|
printf("%10lg", pointlo[k]); |
1711 |
|
|
printf("\n"); |
1712 |
|
|
|
1713 |
|
|
return 0; |
1714 |
|
|
} |
1715 |
|
|
|
1716 |
|
|
#endif |
1717 |
|
|
|
1718 |
|
|
#ifdef TEST_NEXT |
1719 |
|
|
#include <stdio.h> |
1720 |
|
|
|
1721 |
|
|
int main() |
1722 |
|
|
{ |
1723 |
|
|
unsigned i; |
1724 |
|
|
unsigned c1[100], c2[100], pt[100]; |
1725 |
|
|
unsigned nDims, nBytes = 4; |
1726 |
|
|
int stat, findPrev; |
1727 |
|
|
printf("Enter nDims: " ); |
1728 |
|
|
scanf("%u", &nDims); |
1729 |
|
|
|
1730 |
|
|
printf("Enter 1st box corner: "); |
1731 |
|
|
for (i = 0; i < nDims; ++i) |
1732 |
|
|
scanf("%u", &c1[i]); |
1733 |
|
|
printf("Enter 2nd box corner: "); |
1734 |
|
|
for (i = 0; i < nDims; ++i) |
1735 |
|
|
scanf("%u", &c2[i]); |
1736 |
|
|
printf("Enter point: "); |
1737 |
|
|
for (i = 0; i < nDims; ++i) |
1738 |
|
|
scanf("%u", &pt[i]); |
1739 |
|
|
printf("Find prev?: "); |
1740 |
|
|
scanf("%d", &findPrev); |
1741 |
|
|
|
1742 |
|
|
stat = hilbert_nextinbox(nDims, nBytes, 8*nBytes, findPrev, c1, c2, pt); |
1743 |
|
|
|
1744 |
|
|
if (stat) |
1745 |
|
|
for (i = 0; i < nDims; ++i) |
1746 |
|
|
printf("%u ", c1[i]); |
1747 |
|
|
else |
1748 |
|
|
printf("No such point"); |
1749 |
|
|
|
1750 |
|
|
printf("\n"); |
1751 |
|
|
return 0; |
1752 |
|
|
} |
1753 |
|
|
#endif |