1 files changed, 109 insertions, 0 deletions
diff --git a/peanosfc.scm b/peanosfc.scm
new file mode 100644
index 0000000..4a4039a
--- /dev/null
+++ b/peanosfc.scm
@@ -0,0 +1,109 @@
+; "peanospc.scm": Peano space filling mapping
+; Copyright (C) 2005 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)
+
+;;; A. R. Butz.
+;;; Space filling curves and mathematical programming.
+;;; Information and Control, 12:314-330, 1968. 
+
+(define (integer->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))
+	 (set! tets (reverse tets))
+	 (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))))))
+	 tra))))
+
+(define (tet-array->integer tra)
+  (define rank (car (array-dimensions tra)))
+  (define depth (cadr (array-dimensions tra)))
+  (define val 0)
+  (do ((idx 0 (+ 1 idx)))
+      ((>= idx depth) val)
+    (do ((rdx (+ -1 rank) (+ -1 rdx)))
+	((negative? rdx))
+      (set! val (+ (array-ref tra rdx idx) (* 3 val))))))
+
+(define (tet-array->coordinates tra)
+  (define rank (car (array-dimensions tra)))
+  (define depth (cadr (array-dimensions tra)))
+  (do ((rdx (+ -1 rank) (+ -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 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)))
+    (do ((rdx 0 (+ 1 rdx))
+	 (cds coords (cdr cds)))
+	((null? cds))
+      (do ((idx (+ -1 depth) (+ -1 idx))
+	   (scl (car cds) (quotient scl 3)))
+	  ((negative? idx))
+	(array-set! tra (modulo scl 3) rdx idx)))
+    tra))
+
+(define (peano-flip! tra)
+  (define parity 0)
+  (define rank (car (array-dimensions tra)))
+  (define depth (cadr (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)))
+	(array-set! tra (if (odd? tpar) (- 2 v_ij) 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))
+  (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 tra (coordinates->tet-array coords))
+  (peano-flip! tra)
+  (tet-array->integer tra))