1 files changed, 37 insertions, 41 deletions
diff --git a/peanosfc.scm b/peanosfc.scm
index 5cac088..388be9f 100644
--- a/peanosfc.scm
+++ b/peanosfc.scm
@@ -1,5 +1,5 @@
 ; "peanospc.scm": Peano space filling mapping
-; Copyright (C) 2005 Aubrey Jaffer
+; Copyright (C) 2005, 2006 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
@@ -26,24 +26,23 @@
 ;;; Space filling curves and mathematical programming.
 ;;; Information and Control, 12:314-330, 1968. 
 
-(define (natural->tet-array scalar rank)
-  (do ((tets '() (cons (modulo scl 3) tets))
+(define (natural->trit-array scalar rank)
+  (do ((trits '() (cons (modulo scl 3) trits))
        (scl scalar (quotient scl 3)))
       ((zero? scl)
-       (let* ((len (length tets))
-	      (depth (quotient (+ len rank -1) rank)))
+       (let ((depth (quotient (+ (length trits) rank -1) rank)))
 	 (define tra (make-array (A:fixZ8b 0) rank depth))
-	 (set! tets (reverse tets))
+	 (set! trits (reverse trits))
 	 (do ((idx (+ -1 depth) (+ -1 idx)))
 	     ((negative? idx))
 	   (do ((rdx 0 (+ 1 rdx)))
 	       ((>= rdx rank))
-	     (cond ((null? tets))
-		   (else (array-set! tra (car tets) rdx idx)
-			 (set! tets (cdr tets))))))
+	     (cond ((null? trits))
+		   (else (array-set! tra (car trits) rdx idx)
+			 (set! trits (cdr trits))))))
 	 tra))))
 
-(define (tet-array->natural tra)
+(define (trit-array->natural tra)
   (define rank (car (array-dimensions tra)))
   (define depth (cadr (array-dimensions tra)))
   (define val 0)
@@ -53,22 +52,20 @@
 	((negative? rdx))
       (set! val (+ (array-ref tra rdx idx) (* 3 val))))))
 
-(define (tet-array->coordinates tra)
-  (define rank (car (array-dimensions tra)))
+(define (trit-array->natural-coordinates tra)
   (define depth (cadr (array-dimensions tra)))
-  (do ((rdx (+ -1 rank) (+ -1 rdx))
+  (do ((rdx (+ -1 (car (array-dimensions tra))) (+ -1 rdx))
        (lst '() (cons (do ((idx 0 (+ 1 idx))
 			   (val 0 (+ (array-ref tra rdx idx) (* 3 val))))
 			  ((>= idx depth) val))
 		      lst)))
       ((negative? rdx) lst)))
 
-(define (coordinates->tet-array coords)
+(define (natural-coordinates->trit-array coords)
   (define depth (do ((scl (apply max coords) (quotient scl 3))
 		     (dpt 0 (+ 1 dpt)))
 		    ((zero? scl) dpt)))
-  (define rank (length coords))
-  (let ((tra (make-array (A:fixN8b 0) rank depth)))
+  (let ((tra (make-array (A:fixN8b 0) (length coords) depth)))
     (do ((rdx 0 (+ 1 rdx))
 	 (cds coords (cdr cds)))
 	((null? cds))
@@ -82,63 +79,62 @@
   (define parity 0)
   (define rank (car (array-dimensions tra)))
   (define depth (cadr (array-dimensions tra)))
+  (define rra (make-array (A:fixN8b 0) (car (array-dimensions tra))))
   (do ((idx 0 (+ 1 idx)))
       ((>= idx depth))
     (do ((rdx (+ -1 rank) (+ -1 rdx)))
 	((negative? rdx))
-      (let ((v_ij (array-ref tra rdx idx))
-	    (tpar parity))
-	(do ((idx (+ -1 idx) (+ -1 idx)))
-	    ((negative? idx))
-	  (set! tpar (+ (array-ref tra rdx idx) tpar)))
-	(if (odd? tpar) (array-set! tra (- 2 v_ij) rdx idx))
-	(set! parity (modulo (+ parity v_ij) 2))))))
+      (let ((v_ij (array-ref tra rdx idx)))
+	(if (odd? (+ parity (array-ref rra rdx)))
+	    (array-set! tra (- 2 v_ij) rdx idx))
+	(set! parity (modulo (+ v_ij parity) 2))
+	(array-set! rra (modulo (+ v_ij (array-ref rra rdx)) 2) rdx)))))
 
 ;;@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 (natural->peano-coordinates scalar rank)
-  (define tra (natural->tet-array scalar rank))
+  (define tra (natural->trit-array scalar rank))
   (peano-flip! tra)
-  (tet-array->coordinates tra))
+  (trit-array->natural-coordinates tra))
 
 ;;@body
 ;;Returns an exact nonnegative integer corresponding to @1, a list of
 ;;nonnegative integer coordinates.
 (define (peano-coordinates->natural coords)
-  (define tra (coordinates->tet-array coords))
+  (define tra (natural-coordinates->trit-array coords))
   (peano-flip! tra)
-  (tet-array->natural tra))
+  (trit-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)))
+  (define nine^rank (expt 9 rank))
+  (do ((edx 1 (* edx nine^rank))
+       (cdx 1 (* cdx 9)))
       ((>= (quotient edx 2) (abs scalar))
-       (let ((tra (natural->tet-array (+ scalar (quotient edx 2)) rank))
-	     (offset (quotient (expt 3 m) 2)))
+       (let ((tra (natural->trit-array (+ scalar (quotient edx 2)) rank))
+	     (offset (quotient cdx 2)))
 	 (peano-flip! tra)
-	 (map (lambda (k) (* (if (odd? m) -1 1) (- k offset)))
-	      (tet-array->coordinates tra))))))
+	 (map (lambda (k) (- k offset))
+	      (trit-array->natural-coordinates tra))))))
 
 ;;@body
 ;;Returns an exact integer corresponding to @1, a list of integer
 ;;coordinates.
 (define (peano-coordinates->integer coords)
+  (define nine^rank (expt 9 (length coords)))
   (define cobs (apply max (map abs coords)))
-  (let loop ((xpo 1))
-    (define offset (quotient (expt 3 xpo) 2))
+  (let loop ((edx 1) (cdx 1))
+    (define offset (quotient cdx 2))
     (if (>= offset cobs)
-	(let ((tra (coordinates->tet-array
+	(let ((tra (natural-coordinates->trit-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)))))
+	  (- (trit-array->natural tra)
+	     (quotient edx 2)))
+	(loop (* nine^rank edx) (* 9 cdx)))))