Import Upstream version 3a2upstream/3a2

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)))