New upstream version 3b5upstream/3b5upstream

1 files changed, 317 insertions, 0 deletions

diff --git a/rmdsff.scm b/rmdsff.scm
new file mode 100644
index 0000000..b9ee7eb
--- /dev/null
+++ b/rmdsff.scm
@@ -0,0 +1,317 @@
+;;;; "rmdsff.scm" Space-filling functions and their inverses.
+;;; Copyright (C) 2013, 2014 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 'array)
+
+;;@code{(require 'space-filling)}
+;;@ftindex space-filling
+
+;;@ The algorithms and cell properties are described in
+;;@url{http://people.csail.mit.edu/jaffer/Geometry/RMDSFF.pdf}
+
+;;; A cell is an object encapsulating the information about a
+;;; Hamiltonian path on a rectangular grid of side^rank nodes.
+;;; Here are the accessors for a cell:
+(define cell-type caar)
+(define cell-side cadar)
+(define cell-rank caddar)
+(define (cell-index cell crds) (apply array-ref (cadadr cell) crds))
+(define (cell-coords cell t) (vector-ref (caadr cell) t))
+(define (cell-entry cell t) (vector-ref (caddr cell) t))
+(define (cell-exit cell t) (vector-ref (cadddr cell) t))
+(define (cell-rotation cell t) (vector-ref (cadddr (cdr cell)) t))
+
+;;@args type rank side precession
+;;@args type rank side
+;;@args type rank
+;;
+;;@1 must be the symbol @code{diagonal}, @code{adjacent}, or
+;;@code{centered}.  @2 must be an integer larger than 1.  @3, if
+;;present, must be an even integer larger than 1 if @1 is
+;;@code{adjacent} or an odd integer larger than 2 otherwise; @3
+;;defaults to the smallest value.  @4, if present, must be an integer
+;;between 0 and @3^@2-1; it is relevant only when @1 is
+;;@code{diagonal} or @code{centered}.
+;;
+;;@args Hamiltonian-path-vector precession
+;;@args Hamiltonian-path-vector
+;;
+;;@1 must be a vector of @var{side}^@var{rank} lists of @var{rank} of
+;;integers encoding the coordinate positions of a Hamiltonian path on
+;;the @var{rank}-dimensional grid of points starting and ending on
+;;corners of the grid.  The starting corner must be the origin
+;;(all-zero coordinates). If the side-length is even, then the ending
+;;corner must be non-zero in only one coordinate; otherwise, the
+;;ending corner must be the furthest diagonally opposite corner from
+;;the origin.
+;;
+;;@code{make-cell} returns a data object suitable for passing as the
+;;first argument to @code{integer->coordinates} or
+;;@code{coordinates->integer}.
+(define (make-cell arg1 . args)
+  (define (make-serpentine-path rank s)
+    (let loop ((path '(())) (rnk (+ -1 rank)))
+      (if (negative? rnk) path
+	  (loop (let iloop ((seq '()) (sc (+ -1 s)))
+		  (if (negative? sc)
+		      seq
+		      (iloop (append (map (lambda (coords) (cons sc coords))
+					  (if (odd? sc) (reverse path) path))
+				     seq)
+			     (+ -1 sc))))
+		(+ -1 rnk)))))
+  (if (list? arg1) (set! arg1 (list->vector arg1)))
+  (cond
+   ((> (length args) 3)
+    (slib:error 'make-cell 'extra 'arguments 'not 'handled args))
+   ((vector? arg1)
+    (let ((path arg1)
+	  (precession (and (not (null? args)) (car args))))
+      (define frst (vector-ref path 0))
+      (define len (vector-length path))
+      (define s-1 (apply max (apply append (vector->list path))))
+      (let* ((len-1 (+ -1 len))
+	     (last (vector-ref path len-1))
+	     (d (length frst)))
+	;; returns index of first non-zero in LST
+	(define (first-non-zero lst)
+	  (define pos (lambda (n lst)
+			(cond ((zero? (car lst)) (pos (+ 1 n) (cdr lst)))
+			      (else n))))
+	  (pos 0 lst))
+	;; returns the traversal direction of the sub-path.
+	(define (U_e N t)
+	  (if (= t len-1)
+	      (map - (vector-ref path t) (vector-ref path (- t 1)))
+	      (let ((dH_i+1 (map - (vector-ref path (+ t 1)) (vector-ref path t))))
+		(define dotpr (apply + (map * N dH_i+1)))
+		(define csum (apply + dH_i+1))
+		(if (or (and (zero? dotpr) (= 1 csum)) (= dotpr csum -1))
+		    dH_i+1
+		    (map - (vector-ref path t) (vector-ref path (- t 1)))
+		    ))))
+	(define (flip-direction dir cdir)
+	  (map (lambda (px c) (modulo (+ px c) 2)) dir cdir))
+	(define (path-diag? path)
+	  (define prev frst)
+	  (for-each (lambda (lst)
+		      (if (not (= d (length lst)))
+			  (slib:error 'non-uniform 'ranks frst lst)))
+		    (vector->list path))
+	  (for-each (lambda (cs)
+		      (if (not (= 1 (apply + (map abs (map - prev cs)))))
+			  (slib:error 'bad 'step prev cs))
+		      (set! prev cs))
+		    (cdr (vector->list path)))
+	  (cond ((not (zero? (apply + frst))) (slib:error 'strange 'start frst))
+		((not (= d (length last))) (slib:error 'non-uniform 'lengths path))
+		((apply = s-1 last) #t)
+		((and (= s-1 (apply + last))) #f)
+		(else (slib:error 'strange 'net-travel frst last))))
+	(define diag? (path-diag? path))
+	(define entries (make-vector len (vector-ref path 0)))
+	(define exits (make-vector
+		       len (if diag?
+			       (map (lambda (c) (quotient c s-1)) last)
+			       (vector-ref path 1))))
+	(define rotations (make-vector len 0))
+	(define ipath (apply make-array
+			     (A:fixZ32b -1)
+			     (vector->list (make-vector d (+ 1 s-1)))))
+	(define ord 0)
+	(for-each (lambda (coords)
+		    (apply array-set! ipath ord coords)
+		    (set! ord (+ 1 ord)))
+		  (vector->list path))
+	(let lp ((t 1)
+		 (prev-X (if diag?
+			     (map (lambda (c) 1) (vector-ref entries 0))
+			     (vector-ref path 1))))
+	  (cond ((> t len-1)
+		 (if (not (equal? (vector-ref exits len-1)
+				  (map (lambda (c) (quotient c s-1)) last)))
+		     (slib:warn (list (if diag? 'diagonal 'adjacent)
+				      (+ 1 s-1) d precession)
+				'bad 'last 'exit
+				(vector-ref exits len-1) 'should 'be
+				(map (lambda (c) (quotient c s-1)) last)))
+		 (let ((ord 0)
+		       (h (first-non-zero last)))
+		   (for-each
+		    (lambda (coords)
+		      (vector-set! rotations
+				   ord
+				   (modulo
+				    (cond ((not diag?)
+					   (- h (first-non-zero
+						 (map -
+						      (vector-ref exits ord)
+						      (vector-ref entries ord)))))
+					  (precession (+ precession ord))
+					  (else 0))
+				    d))
+		      (set! ord (+ 1 ord)))
+		    (vector->list path)))
+		 (list (list (if diag? 'diagonal 'adjacent)
+			     (+ 1 s-1)
+			     d
+			     precession)
+		       (list path ipath)
+		       entries
+		       exits
+		       rotations
+		       ))
+		(else
+		 (let ((N (flip-direction
+			   prev-X
+			   (map - (vector-ref path t)
+				(vector-ref path (- t 1))))))
+		   (define X (if diag?
+				 (map (lambda (tn) (- 1 tn)) N)
+				 (flip-direction N (U_e N t))))
+		   (vector-set! entries t N)
+		   (vector-set! exits t X)
+		   (lp (+ 1 t) X))))))))
+   ((< (car args) 2)
+    (slib:error 'make-cell 'rank 'too 'small (car args)))
+   (else
+    (case arg1
+      ((center centered)
+       (let ((cell (make-cell (make-serpentine-path
+			       (car args)
+			       (if (null? (cdr args)) 3 (cadr args)))
+			      (if (= 3 (length args)) (caddr args) #f))))
+	 (if (not (eq? 'diagonal (cell-type cell)))
+	     (slib:error 'make-cell 'centered 'must 'be 'diagonal (car cell)))
+	 (set-car! (car cell) 'centered)
+	 cell))
+      ((diagonal opposite)
+       (make-cell (make-serpentine-path
+		   (car args)
+		   (if (null? (cdr args)) 3 (cadr args)))
+		  (if (= 3 (length args)) (caddr args) #f)))
+      ((adjacent)
+       (make-cell (make-serpentine-path
+		   (car args)
+		   (if (null? (cdr args)) 2 (cadr args)))
+		  (if (= 3 (length args)) (caddr args) #f)))
+      (else
+       (slib:error 'make-cell 'unknown 'cell 'type arg1))))))
+
+;;@ Hilbert, Peano, and centered Peano cells are generated
+;;respectively by:
+;;@example
+;;(make-cell 'adjacent @var{rank} 2)   ; Hilbert
+;;(make-cell 'diagonal @var{rank} 3)   ; Peano
+;;(make-cell 'centered @var{rank} 3)   ; centered Peano
+;;@end example
+
+;; Positive k rotates left
+(define (rotate-list lst k)
+  (define len (length lst))
+  (cond ((<= len 1) lst)
+	(else
+	 (set! k (modulo k len))
+	 (if (zero? k)
+	     lst
+	     (let ((ans (list (car lst))))
+	       (do ((ls (cdr lst) (cdr ls))
+		    (tail ans (cdr tail))
+		    (k (+ -1 k) (+ -1 k)))
+		   ((<= k 0)
+		    (append ls ans))
+		 (set-cdr! tail (list (car ls)))))))))
+
+;;@ In the conversion procedures, if the cell is @code{diagonal} or
+;;@code{adjacent}, then the coordinates and scalar must be nonnegative
+;;integers.  If @code{centered}, then the integers can be negative.
+
+;;@body
+;;@0 converts the integer @2 to a list of coordinates according to @1.
+(define (integer->coordinates cell u)
+  (define umag (case (cell-type cell)
+		 ((centered) (* (abs u) 2))
+		 (else u)))
+  (define d (cell-rank cell))
+  (define s (cell-side cell))
+  (let* ((s^d (expt s d))
+	 (s^d^2 (expt s^d d)))
+    (define (align V t sde)
+      (map (lambda (Vj Nj) (if (zero? Nj) Vj (- sde Vj)))
+	   (rotate-list V (cell-rotation cell t))
+	   (cell-entry cell t)))
+    (define (rec u m w)
+      (if (positive? w)
+	  (let ((t (quotient u m)))
+	    (map +
+		 (map (lambda (y) (* y w)) (cell-coords cell t))
+		 (align (rec (modulo u m) (quotient m s^d) (quotient w s))
+			t
+			(- w 1))))
+	  (cell-coords cell 0)))
+    (do ((uscl 1 (* uscl s^d^2))
+	 (cscl 1 (* cscl s^d)))
+	((> uscl umag)
+	 (case (cell-type cell)
+	   ((centered)
+	    (let ((cscl/2 (quotient cscl 2)))
+	      (map (lambda (c) (- c cscl/2))
+		   (rec (+ u (quotient uscl 2))
+			(quotient uscl s^d)
+			(quotient cscl s)))))
+	   (else (rec u (quotient uscl s^d) (quotient cscl s))))))))
+
+;;@body
+;;@0 converts the list of coordinates @2 to an integer according to @1.
+(define (coordinates->integer cell V)
+  (define maxc (case (cell-type cell)
+		 ((centered) (* 2 (apply max (map abs V))))
+		 (else (apply max V))))
+  (define d (cell-rank cell))
+  (define s (cell-side cell))
+  (let* ((s^d (expt s d))
+	 (s^d^2 (expt s^d d)))
+    (define (align^-1 V t sde)
+      (rotate-list (map (lambda (Vj Nj) (if (zero? Nj) Vj (- sde Vj)))
+			V
+			(cell-entry cell t))
+		   (- (cell-rotation cell t))))
+    (define (rec u V w)
+      (if (positive? w)
+	  (let ((dig (cell-index cell (map (lambda (c) (quotient c w)) V))))
+	    (rec (+ dig (* s^d u))
+		 (align^-1 (map (lambda (cx) (modulo cx w)) V)
+			   dig
+			   (- w 1))
+		 (quotient w s)))
+	  u))
+    (do ((uscl 1 (* uscl s^d^2))
+	 (cscl 1 (* cscl s^d)))
+	((> cscl maxc)
+	 (case (cell-type cell)
+	   ((centered)
+	    (let ((cscl/2 (quotient cscl 2)))
+	      (- (rec 0
+		      (map (lambda (c) (+ c cscl/2)) V)
+		      (quotient cscl s))
+		 (quotient uscl 2))))
+	   (else (rec 0 V (quotient cscl s))))))))
+
+;;@var{coordinates->integer} and @var{integer->coordinates} are
+;;inverse functions when passed the same @var{cell} argument.