; "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
;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 Hilbert
;;@cindex Space-Filling
;;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
;;arbitrarily large @var{n}-dimensional cube with its corner at the
;;origin and all coordinates are non-negative.
;;
;;@noindent
;;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



;;@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.
;;
;;@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
    (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))))))

;;@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 . nbits)
  (define rank (length coords))
  (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)
		  (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)))