diff options
Diffstat (limited to 'phil-spc.scm')
-rw-r--r-- | phil-spc.scm | 242 |
1 files changed, 197 insertions, 45 deletions
diff --git a/phil-spc.scm b/phil-spc.scm index 3372ce6..65863da 100644 --- a/phil-spc.scm +++ b/phil-spc.scm @@ -1,5 +1,5 @@ -; "phil-spc.scm": Peano-Hilbert space filling mapping -; Copyright (c) 2003 Aubrey Jaffer +; "phil-spc.scm": 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 @@ -23,11 +23,12 @@ ;;@ftindex hilbert-fill ;; ;;@noindent -;;@cindex Peano ;;@cindex Hilbert ;;@cindex Space-Filling -;;The @dfn{Peano-Hilbert Space-Filling Curve} is a one-to-one mapping +;;The @dfn{Hilbert Space-Filling Curve} is a one-to-one mapping ;;between a unit line segment and an @var{n}-dimensional unit cube. +;;This implementation treats the nonnegative integers either as +;;fractional bits of a given width or as nonnegative integers. ;; ;;@noindent ;;The integer procedures map the non-negative integers to an @@ -35,60 +36,211 @@ ;;origin and all coordinates are non-negative. ;; ;;@noindent -;;For any exact nonnegative integers @var{scalar} and @var{rank}, +;;For any exact nonnegative integer @var{scalar} and exact integer +;;@var{rank} > 2, ;; ;;@example ;;(= @var{scalar} (hilbert-coordinates->integer ;; (integer->hilbert-coordinates @var{scalar} @var{rank}))) ;; @result{} #t ;;@end example +;; +;;When treating integers as @var{k} fractional bits, +;; +;;@example +;;(= @var{scalar} (hilbert-coordinates->integer +;; (integer->hilbert-coordinates @var{scalar} @var{rank} @var{k})) @var{k}) +;; @result{} #t +;;@end example -;;@body + + +;;@args scalar rank ;;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)) +;; +;;@args scalar rank k +;; +;;@1 must be a nonnegative integer of no more than +;;@code{@2*@var{k}} bits. +;; +;;@0 Returns a list of @2 @var{k}-bit nonnegative integer +;;coordinates corresponding to exact non-negative integer @1. The +;;curves generated by @0 have the same alignment independent of +;;@var{k}. +(define (integer->hilbert-coordinates scalar rank . nbits) + (define igry (integer->gray-code scalar)) + (define rnkmsk (lognot (ash -1 rank))) + (define rnkhib (ash 1 (+ -1 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)))))) + (if (null? nbits) + (let ((rank^2 (* rank rank))) + (* (quotient (+ -1 rank^2 (integer-length scalar)) rank^2) + rank^2)) + (* rank (car nbits)))) + (do ((bdxn (- rank rank*nbits) (+ rank bdxn)) + (chnk (ash igry (- rank rank*nbits)) + (logxor rnkhib (logand (ash igry (+ rank bdxn)) rnkmsk))) + (rotation 0 (modulo (+ (log2-binary-factors chnk) 2 rotation) rank)) + (flipbit 0 (ash 1 rotation)) + (lst '() (cons (logxor flipbit (rotate-bit-field chnk rotation 0 rank)) + lst))) + ((positive? bdxn) + (map gray-code->integer (delaminate-list rank (reverse lst)))))) -;;@body +;;@args coords +;;@args coords k ;;Returns an exact non-negative integer corresponding to @1, a list ;;of non-negative integer coordinates. -(define (hilbert-coordinates->integer coords) +(define (hilbert-coordinates->integer coords . nbits) (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) + (let ((lst (delaminate-list rank (map integer->gray-code coords))) + (rnkhib (ash 1 (+ -1 rank)))) + (define (loop lst rotation flipbit scalar) + (if (null? lst) + (gray-code->integer scalar) + (let ((chnk (rotate-bit-field (logxor flipbit (car lst)) + (- rotation) 0 rank))) + (loop (cdr lst) + (modulo (+ (log2-binary-factors chnk) 2 rotation) rank) (ash 1 rotation) - (+ chnk (ash scalar rank)))))) - (loop (- rank rank*nbits) 0 0 0))) + (logior (logxor rnkhib chnk) (ash scalar rank)))))) + (loop (cdr lst) + (modulo (+ (log2-binary-factors (car lst)) 2) rank) + 1 + (car lst)))) + +;;@subsubsection Gray code +;; +;;@cindex Gray code +;;@noindent +;;A @dfn{Gray code} is an ordering of non-negative integers in which +;;exactly one bit differs between each pair of successive elements. There +;;are multiple Gray codings. An n-bit Gray code corresponds to a +;;Hamiltonian cycle on an n-dimensional hypercube. +;; +;;@noindent +;;Gray codes find use communicating incrementally changing values between +;;asynchronous agents. De-laminated Gray codes comprise the coordinates +;;of Hilbert space-filling curves. +;; +;; +;;@defun integer->gray-code k +;;Converts @var{k} to a Gray code of the same @code{integer-length} as +;;@var{k}. +;; +;;@defunx gray-code->integer k +;;Converts the Gray code @var{k} to an integer of the same +;;@code{integer-length} as @var{k}. +;; +;;For any non-negative integer @var{k}, +;;@example +;;(eqv? k (gray-code->integer (integer->gray-code k))) +;;@end example +;;@end defun +(define (integer->gray-code k) + (logxor k (arithmetic-shift k -1))) +(define (gray-code->integer k) + (if (negative? k) + (slib:error 'gray-code->integer 'negative? k) + (let ((kln (integer-length k))) + (do ((d 1 (* d 2)) + (ans (logxor k (arithmetic-shift k -1)) ; == (integer->gray-code k) + (logxor ans (arithmetic-shift ans (* d -2))))) + ((>= (* 2 d) kln) ans))))) + +(define (grayter k1 k2) + (define kl1 (integer-length k1)) + (define kl2 (integer-length k2)) + (if (eqv? kl1 kl2) + (> (gray-code->integer k1) (gray-code->integer k2)) + (> kl1 kl2))) + +;;@defun = k1 k2 +;;@defunx gray-code<? k1 k2 +;;@defunx gray-code>? k1 k2 +;;@defunx gray-code<=? k1 k2 +;;@defunx gray-code>=? k1 k2 +;;These procedures return #t if their Gray code arguments are +;;(respectively): equal, monotonically increasing, monotonically +;;decreasing, monotonically nondecreasing, or monotonically nonincreasing. +;; +;;For any non-negative integers @var{k1} and @var{k2}, the Gray code +;;predicate of @code{(integer->gray-code k1)} and +;;@code{(integer->gray-code k2)} will return the same value as the +;;corresponding predicate of @var{k1} and @var{k2}. +;;@end defun +(define (gray-code<? k1 k2) + (not (or (eqv? k1 k2) (grayter k1 k2)))) +(define (gray-code<=? k1 k2) + (or (eqv? k1 k2) (not (grayter k1 k2)))) +(define (gray-code>? k1 k2) + (and (not (eqv? k1 k2)) (grayter k1 k2))) +(define (gray-code>=? k1 k2) + (or (eqv? k1 k2) (grayter k1 k2))) + +;;@subsubsection Bitwise Lamination +;;@cindex lamination + +;;@args k1 @dots{} +;;Returns an integer composed of the bits of @var{k1} @dots{} interlaced +;;in argument order. Given @var{k1}, @dots{} @var{kn}, the n low-order +;;bits of the returned value will be the lowest-order bit of each +;;argument. +;; +;;@args count k +;;Returns a list of @var{count} integers comprised of every @var{count}h +;;bit of the integer @var{k}. +;; +;;@example +;;(map (lambda (k) (number->string k 2)) +;; (bitwise-delaminate 4 #x7654)) +;; @result{} ("0" "1111" "1100" "1010") +;;(number->string (bitwise-laminate 0 #b1111 #b1100 #b1010) 16) +;; @result{} "7654" +;@end example +;; +;;For any non-negative integers @var{k} and @var{count}: +;;@example +;;(eqv? k (bitwise-laminate (bitwise-delaminate count k))) +;;@end example +(define (bitwise-laminate . ks) + (define nks (length ks)) + (define nbs (apply max (map integer-length ks))) + (do ((kdx (+ -1 nbs) (+ -1 kdx)) + (ibs 0 (+ (list->integer (map (lambda (k) (logbit? kdx k)) ks)) + (arithmetic-shift ibs nks)))) + ((negative? kdx) ibs))) +(define (bitwise-delaminate count k) + (define nbs (* count (+ 1 (quotient (integer-length k) count)))) + (do ((kdx (- nbs count) (- kdx count)) + (lst (vector->list (make-vector count 0)) + (map (lambda (k bool) (+ (if bool 1 0) (arithmetic-shift k 1))) + lst + (integer->list (arithmetic-shift k (- kdx)) count)))) + ((negative? kdx) lst))) + +;;@body +;; +;;Returns a list of @var{count} integers comprised of the @var{j}th +;;bit of the integers @var{ks} where @var{j} ranges from @var{count}-1 +;;to 0. +;; +;;@example +;;(map (lambda (k) (number->string k 2)) +;; (delaminate-list 4 '(7 6 5 4 0 0 0 0))) +;; @result{} ("0" "11110000" "11000000" "10100000") +;;@end example +;; +;;@0 is its own inverse: +;;@example +;;(delaminate-list 8 (delaminate-list 4 '(7 6 5 4 0 0 0 0))) +;; @result{} (7 6 5 4 0 0 0 0) +;;@end example +(define (delaminate-list count ks) + (define nks (length ks)) + (do ((kdx 0 (+ 1 kdx)) + (lst '() (cons (list->integer (map (lambda (k) (logbit? kdx k)) ks)) + lst))) + ((>= kdx count) lst))) |