; "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)

;;@code{(require 'peano-fill)}
;;@ftindex peano-fill

;;; A. R. Butz.
;;; 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))
       (scl scalar (quotient scl 3)))
      ((zero? scl)
       (let* ((len (length tets))
	      (depth (quotient (+ len rank -1) rank)))
	 (define tra (make-array (A:fixZ8b 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->natural 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 (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->natural coords)
  (define tra (coordinates->tet-array coords))
  (peano-flip! 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)))))