Import Upstream version 3a3upstream/3a3

1 files changed, 43 insertions, 8 deletions

diff --git a/peanosfc.scm b/peanosfc.scm
index 4a4039a..5cac088 100644
--- a/peanosfc.scm
+++ b/peanosfc.scm
@@ -19,17 +19,20 @@
 
 (require 'array)
 
+;;@code{(require 'peano-fill)}
+;;@ftindex peano-fill
+
 ;;; A. R. Butz.
 ;;; Space filling curves and mathematical programming.
 ;;; Information and Control, 12:314-330, 1968. 
 
-(define (integer->tet-array scalar rank)
+(define (natural->tet-array scalar rank)
   (do ((tets '() (cons (modulo scl 3) tets))
        (scl scalar (quotient scl 3)))
       ((zero? scl)
        (let* ((len (length tets))
 	      (depth (quotient (+ len rank -1) rank)))
-	 (define tra (make-array (A:fixN8b 0) rank depth))
+	 (define tra (make-array (A:fixZ8b 0) rank depth))
 	 (set! tets (reverse tets))
 	 (do ((idx (+ -1 depth) (+ -1 idx)))
 	     ((negative? idx))
@@ -40,7 +43,7 @@
 			 (set! tets (cdr tets))))))
 	 tra))))
 
-(define (tet-array->integer tra)
+(define (tet-array->natural tra)
   (define rank (car (array-dimensions tra)))
   (define depth (cadr (array-dimensions tra)))
   (define val 0)
@@ -88,22 +91,54 @@
 	(do ((idx (+ -1 idx) (+ -1 idx)))
 	    ((negative? idx))
 	  (set! tpar (+ (array-ref tra rdx idx) tpar)))
-	(array-set! tra (if (odd? tpar) (- 2 v_ij) v_ij) rdx idx)
+	(if (odd? tpar) (array-set! tra (- 2 v_ij) rdx idx))
 	(set! parity (modulo (+ parity v_ij) 2))))))
 
 ;;@body
 ;;Returns a list of @2 nonnegative integer coordinates corresponding
 ;;to exact nonnegative integer @1.  The lists returned by @0 for @1
 ;;arguments 0 and 1 will differ in the first element.
-(define (integer->peano-coordinates scalar rank)
-  (define tra (integer->tet-array scalar rank))
+(define (natural->peano-coordinates scalar rank)
+  (define tra (natural->tet-array scalar rank))
   (peano-flip! tra)
   (tet-array->coordinates tra))
 
 ;;@body
 ;;Returns an exact nonnegative integer corresponding to @1, a list of
 ;;nonnegative integer coordinates.
-(define (peano-coordinates->integer coords)
+(define (peano-coordinates->natural coords)
   (define tra (coordinates->tet-array coords))
   (peano-flip! tra)
-  (tet-array->integer tra))
+  (tet-array->natural tra))
+
+;;@body
+;;Returns a list of @2 integer coordinates corresponding to exact
+;;integer @1.  The lists returned by @0 for @1 arguments 0 and 1 will
+;;differ in the first element.
+(define (integer->peano-coordinates scalar rank)
+  (define three^rank (expt 3 rank))
+  (do ((edx 1 (* edx three^rank))
+       (m 0 (+ 1 m)))
+      ((>= (quotient edx 2) (abs scalar))
+       (let ((tra (natural->tet-array (+ scalar (quotient edx 2)) rank))
+	     (offset (quotient (expt 3 m) 2)))
+	 (peano-flip! tra)
+	 (map (lambda (k) (* (if (odd? m) -1 1) (- k offset)))
+	      (tet-array->coordinates tra))))))
+
+;;@body
+;;Returns an exact integer corresponding to @1, a list of integer
+;;coordinates.
+(define (peano-coordinates->integer coords)
+  (define cobs (apply max (map abs coords)))
+  (let loop ((xpo 1))
+    (define offset (quotient (expt 3 xpo) 2))
+    (if (>= offset cobs)
+	(let ((tra (coordinates->tet-array
+		    (map (lambda (elt) (+ elt offset))
+			 coords))))
+	  (peano-flip! tra)
+	  ((if (odd? xpo) - +)
+	   (- (tet-array->natural tra)
+	      (quotient (expt 3 (* (length coords) xpo)) 2))))
+	(loop (+ 1 xpo)))))