1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
; "phil-spc.scm": Peano-Hilbert space filling mapping
; Copyright (c) 2003 Aubrey Jaffer
;
;Permission to copy this software, to modify it, to redistribute it,
;to distribute modified versions, and to use it for any purpose is
;granted, subject to the following restrictions and understandings.
;
;1. Any copy made of this software must include this copyright notice
;in full.
;
;2. I have made no warranty or representation that the operation of
;this software will be error-free, and I am under no obligation to
;provide any services, by way of maintenance, update, or otherwise.
;
;3. In conjunction with products arising from the use of this
;material, there shall be no use of my name in any advertising,
;promotional, or sales literature without prior written consent in
;each case.
(require 'logical)
;;@code{(require 'hilbert-fill)}
;;@ftindex hilbert-fill
;;
;;@noindent
;;@cindex Peano
;;@cindex Hilbert
;;@cindex Space-Filling
;;The @dfn{Peano-Hilbert Space-Filling Curve} is a one-to-one mapping
;;between a unit line segment and an @var{n}-dimensional unit cube.
;;
;;@noindent
;;The integer procedures map the non-negative integers to an
;;arbitrarily large @var{n}-dimensional cube with its corner at the
;;origin and all coordinates are non-negative.
;;
;;@noindent
;;For any exact nonnegative integers @var{scalar} and @var{rank},
;;
;;@example
;;(= @var{scalar} (hilbert-coordinates->integer
;; (integer->hilbert-coordinates @var{scalar} @var{rank})))
;; @result{} #t
;;@end example
;;@body
;;Returns a list of @2 integer coordinates corresponding to exact
;;non-negative integer @1. The lists returned by @0 for @1 arguments
;;0 and 1 will differ in the first element.
(define (integer->hilbert-coordinates scalar rank)
(define ndones (logical:ones rank))
(define rank*nbits
(let ((rank^2 (* rank rank)))
(* (quotient (+ -1 rank^2 (integer-length scalar)) rank^2)
rank^2)))
(let ((nthbits (quotient (logical:ones rank*nbits) ndones)))
(define igry (logxor (integer->gray-code scalar) (ash nthbits -1)))
(do ((bdxn (- rank rank*nbits) (+ rank bdxn))
(chnk (logand (ash igry (- rank rank*nbits)) ndones)
(logand (ash igry (+ rank bdxn)) ndones))
(rotation 0 (modulo (+ (integer-length (logand (- chnk) chnk))
1 rotation)
rank))
(flipbit 0 (ash 1 rotation))
(bignum 0 (+ (logxor flipbit (logical:rotate chnk rotation rank))
(ash bignum rank))))
((positive? bdxn)
(map gray-code->integer (bitwise:delaminate rank bignum))))))
;;@body
;;Returns an exact non-negative integer corresponding to @1, a list
;;of non-negative integer coordinates.
(define (hilbert-coordinates->integer coords)
(define rank (length coords))
(define bignum (apply bitwise:laminate (map integer->gray-code coords)))
(let ((rank*nbits
(* (quotient (+ -1 rank (integer-length (apply max coords))) rank)
rank rank))
(ndones (logical:ones rank)))
(define nthbits (quotient (logical:ones rank*nbits) ndones))
(define (loop bdxn rotation flipbit scalar)
(if (positive? bdxn)
(gray-code->integer (logxor scalar (ash nthbits -1)))
(let ((chnk (logical:rotate
(logxor flipbit (logand ndones (ash bignum bdxn)))
(- rotation)
rank)))
(loop (+ rank bdxn)
(modulo (+ (integer-length (logand (- chnk) chnk))
1 rotation)
rank)
(ash 1 rotation)
(+ chnk (ash scalar rank))))))
(loop (- rank rank*nbits) 0 0 0)))
|
