(declare (fixnum))

;;; The following macro is used to determine whether certain keyword arguments

;;; were omitted.  It is specific to Gambit-C's compiler; note that it's not defined

;;; in the Gambit-C interpreter.  Redefine it for other schemes.

(##define-macro (macro-absent-obj)

  (##namespace ("c#" absent-object))

  `',absent-object)

;;; This is Gambit-specific to make it a bit faster

(define (exact-integer? x)

  (or (##fixnum? x)

      (##bignum? x)))

;;; We need a multi-argument every, but not something as fancy as in Olin Shivers

;;; list library.  (Shiver's version works fine, though, for our purposes.)

(define (every pred list . lists)

  (if (pair? lists)

      (let loop ((lists (cons list lists)))

	(or (null? (car lists))

	    (and (apply pred (map car lists))

		 (loop (map cdr lists)))))

      (let loop ((list list))

	(or (null? list)

	    (and (pred (car list))

		 (loop (cdr list)))))))

;;; An Interval is a cross product of multi-indices

;;; [l_0,u_0) x [l_1,u_1) x ... x [l_n-1,u_n-1)

;;; where l_i < u_i for 0 <= i < n, and n > 0 is the dimension of the Interval

(define-structure Interval

  lower-bounds            ;; a vector of exact integers l_0,...,l_n-1

  upper-bounds)           ;; a vector of exact integers u_0,...,u_n-1

(define (build-Interval lower-bounds upper-bounds)

  (or (and (not (vector? lower-bounds))

	   (error "build-Interval: lower-bounds must be a vector: " lower-bounds))

      (and (not (vector? upper-bounds))

	   (error "build-Interval: upper-bounds must be a vector: " upper-bounds))

      (and (not (= (vector-length lower-bounds) (vector-length upper-bounds)))

	   (error "build-Interval: lower-bounds and upper-bounds must be the same length: " lower-bounds upper-bounds))

      (and (not (< 0 (vector-length lower-bounds)))

	   (error "build-Interval: lower-bounds and upper-bounds must be nonempty vectors: " lower-bounds upper-bounds))

      (let ((lower-bounds-list (vector->list lower-bounds))

	    (upper-bounds-list (vector->list upper-bounds)))

	(or (and (not (every exact-integer? lower-bounds-list))

		 (error "build-Interval: All lower-bounds must be exact integers: " lower-bounds))

	    (and (not (every exact-integer? upper-bounds-list))

		 (error "build-Interval: All upper-bounds must be exact integers: " upper-bounds))

	    (and (not (every (lambda (x y) (< x y)) lower-bounds-list upper-bounds-list))

		 (error "build-Interval: Each lower-bound must be less than the associated upper-bound: " lower-bounds upper-bounds))

	    (make-Interval lower-bounds upper-bounds)))))

(define (##Interval-dimension interval)

  (vector-length (Interval-lower-bounds interval)))

(define (Interval-dimension interval)

  (or (and (not (Interval? interval))

	   (error "Interval-dimension: argument is not an interval: " interval))

      (##Interval-dimension interval)))

(define (##Interval-lower-bound interval i)

  (vector-ref (Interval-lower-bounds interval) i))

(define (##Interval-upper-bound interval i)

  (vector-ref (Interval-upper-bounds interval) i))

(define (Interval-lower-bound interval i)

  (or (and (not (Interval? interval))

	   (error "Interval-lower-bound: argument is not an interval: " interval i))

      (and (not (exact-integer? i))

	   (error "Interval-lower-bound: argument is not a exact integer: " interval i))

      (and (not (< -1 i (##Interval-dimension interval)))

	   (error "Interval-lower-bound: index is not between 0 (inclusive) and (Interval-dimension interval) (exclusive): " interval i))

      (vector-ref (Interval-lower-bounds interval) i)))

(define (Interval-upper-bound interval i)

  (or (and (not (Interval? interval))

	   (error "Interval-upper-bound: argument is not an interval: " interval i))

      (and (not (exact-integer? i))

	   (error "Interval-upper-bound: argument is not a exact integer: " interval i))

      (and (not (< -1 i (##Interval-dimension interval)))

	   (error "Interval-upper-bound: index is not between 0 (inclusive) and (Interval-dimension interval) (exclusive): " interval i))

      (vector-ref (Interval-upper-bounds interval) i)))

(define (##Interval-lower-bounds->list interval)

  (vector->list (Interval-lower-bounds interval)))

(define (Interval-lower-bounds->list interval)

  (or (and (not (Interval? interval))

	   (error "Interval-lower-bounds->list: argument is not an interval: " interval))

      (##Interval-lower-bounds->list interval)))

(define (##Interval-upper-bounds->list interval)

  (vector->list (Interval-upper-bounds interval)))

(define (Interval-upper-bounds->list interval)

  (or (and (not (Interval? interval))

	   (error "Interval-upper-bounds->list: argument is not an interval: " interval))

      (##Interval-upper-bounds->list interval)))

(define (Interval-curry interval left-dimension)

  (declare (fixnum))

  (or (and (not (Interval? interval))

	   (error "Interval-curry argument is not an interval: " interval left-dimension))

      (and (not (exact-integer? left-dimension))

	   (error "Interval-curry argument is not a exact integer: " interval left-dimension))

      (and (not (< 0 left-dimension (##Interval-dimension interval)))

	   (error "Interval-curry: argument is not between 0 and (Interval-dimension interval) (exclusive): " interval left-dimension))

      (let ((n (##Interval-dimension interval)))

	(let ((lower-bounds (Interval-lower-bounds interval))

	      (upper-bounds (Interval-upper-bounds interval))

	      (left-lower-bounds (make-vector left-dimension))

	      (left-upper-bounds (make-vector left-dimension))

	      (right-lower-bounds (make-vector (- n left-dimension)))

	      (right-upper-bounds (make-vector (- n left-dimension))))

	  (do ((i 0 (+ i 1)))

	      ((= i left-dimension)

	       (do ((i i (+ i 1)))

		   ((= i n)

		    (values (make-Interval left-lower-bounds

					   left-upper-bounds)

			    (make-Interval right-lower-bounds

					   right-upper-bounds)))

		 (vector-set! right-lower-bounds (- i left-dimension) (vector-ref lower-bounds i))

		 (vector-set! right-upper-bounds (- i left-dimension) (vector-ref upper-bounds i))))

	    (vector-set! left-lower-bounds i (vector-ref lower-bounds i))

	    (vector-set! left-upper-bounds i (vector-ref upper-bounds i)))))))

(define (Interval-distinguish-one interval index)

    (or (and (not (Interval? interval))

	   (error "Interval-distinguish-one argument is not an interval: " interval index))

      (and (not (exact-integer? index))

	   (error "Interval-distinguish-one argument is not a exact integer: " interval index))

      (and (not (< 1 (##Interval-dimension interval)))

	   (error "Interval-distinguish-one: The dimension of the argument is not greater than one: " interval index))

      (and (not (< -1 index (##Interval-dimension interval)))

	   (error "Interval-distinguish-one: argument is not between 0 (inclusive) and (Interval-dimension interval) (exclusive): " interval index))

      (let ((n (##Interval-dimension interval)))

	(let* ((lower-bounds (Interval-lower-bounds interval))

	       (upper-bounds (Interval-upper-bounds interval))

	       (left-lower-bounds (make-vector (- n 1)))

	       (left-upper-bounds (make-vector (- n 1))))

	  (let loop ((i 0)

		     (j 0))

	    (cond ((= i n)

		   (values (build-Interval left-lower-bounds

					   left-upper-bounds)

			   (build-Interval (vector (vector-ref lower-bounds index))

					   (vector (vector-ref upper-bounds index)))))

		  ((= i index)

		   (loop (+ i 1)

			 j))

		  (else

		   (vector-set! left-lower-bounds j (vector-ref lower-bounds i))

		   (vector-set! left-upper-bounds j (vector-ref upper-bounds i))

		   (loop (+ i 1)

			 (+ j 1)))))))))

(define (Interval-volume interval)

  (or (and (not (Interval? interval))

	   (error "Interval-volume: argument is not an interval: " interval))

      (do ((i (- (##Interval-dimension interval) 1) (- i 1))

	   (result 1 (* result (- (##Interval-upper-bound interval i)

				  (##Interval-lower-bound interval i)))))

	  ((< i 0) result))))

(define (Interval= interval1 interval2)

  (or (and (not (and (Interval? interval1) (Interval? interval2)))

	   (error "Interval=: Not all arguments are intervals: " interval1 interval2))

      (and (= (##Interval-dimension interval1) (##Interval-dimension interval2))

	   (let loop ((i (- (##Interval-dimension interval1) 1)))

	     (or (< i 0)

		 (and (= (##Interval-upper-bound interval1 i)

			 (##Interval-upper-bound interval2 i))

		      (= (##Interval-lower-bound interval1 i)

			 (##Interval-lower-bound interval2 i))

		      (loop (- i 1))))))))

(define (Interval-subset? interval1 interval2)

  (or (and (not (and (Interval? interval1) (Interval? interval2)))

	   (error "Interval-subset?: Not all arguments are intervals: " interval1 interval2))

      (and (= (##Interval-dimension interval1) (##Interval-dimension interval2))

	   (let loop ((i (- (##Interval-dimension interval1) 1)))

	     (or (< i 0)

		 (and (>= (##Interval-lower-bound interval1 i)

			  (##Interval-lower-bound interval2 i))

		      (<= (##Interval-upper-bound interval1 i)

			  (##Interval-upper-bound interval2 i))

		      (loop (- i 1))))))))

(define (Interval-contains-multi-index? interval i

					#!optional

					(j (macro-absent-obj))

					(k (macro-absent-obj))

					(l (macro-absent-obj))

					#!rest

					multi-index-tail)

  (or (and (not (Interval? interval))

	   (error "Interval-contains-multi-index?: argument is not an Interval: " interval))

      (cond ((eq? j (macro-absent-obj))

	     (or (and (not (exact-integer? i))

		      (error "Interval-contains-multi-index?: multi-index component is not an exact integer: " i))

		 (and (not (= (##Interval-dimension interval) 1))

		      (error "Interval-contains-multi-index?: dimension of interval does not match number of arguments: " interval i))

		 (and (<=   (##Interval-lower-bound interval 0) i)

		      (<  i (##Interval-upper-bound interval 0)))))

	    ((eq? k (macro-absent-obj))

	     (or (and (not (and (exact-integer? i)

				(exact-integer? j)))

		      (error "Interval-contains-multi-index?: at least one multi-index component is not an exact integer: " i j))

		 (and (not (= (##Interval-dimension interval) 2))

		      (error "Interval-contains-multi-index?: dimension of interval does not match number of arguments: " interval i j))

		 (and (<=   (##Interval-lower-bound interval 0) i)

		      (<  i (##Interval-upper-bound interval 0))

		      (<=   (##Interval-lower-bound interval 1) j)

		      (<  j (##Interval-upper-bound interval 1)))))

	    ((eq? l (macro-absent-obj))

	     (or (and (not (and (exact-integer? i)

				(exact-integer? j)

				(exact-integer? k)))

		      (error "Interval-contains-multi-index?: at least one multi-index component is not an exact integer: " i j k))

		 (and (not (= (##Interval-dimension interval) 3))

		      (error "Interval-contains-multi-index?: dimension of interval does not match number of arguments: " interval i j k))

		 (and (<=   (##Interval-lower-bound interval 0) i)

		      (<  i (##Interval-upper-bound interval 0))

		      (<=   (##Interval-lower-bound interval 1) j)

		      (<  j (##Interval-upper-bound interval 1))

		      (<=   (##Interval-lower-bound interval 2) k)

		      (<  k (##Interval-upper-bound interval 2)))))

	    ((null? multi-index-tail)

	     (or (and (not (and (exact-integer? i)

				(exact-integer? j)

				(exact-integer? k)

				(exact-integer? l)))

		      (error "Interval-contains-multi-index?: at least one multi-index component is not an exact integer: " i j k l))

		 (and (not (= (##Interval-dimension interval) 4))

		      (error "Interval-contains-multi-index?: dimension of interval does not match number of arguments: " interval i j k l))

		 (and (<=   (##Interval-lower-bound interval 0) i)

		      (<  i (##Interval-upper-bound interval 0))

		      (<=   (##Interval-lower-bound interval 1) j)

		      (<  j (##Interval-upper-bound interval 1))

		      (<=   (##Interval-lower-bound interval 2) k)

		      (<  k (##Interval-upper-bound interval 2))

		      (<=   (##Interval-lower-bound interval 3) l)

		      (<  l (##Interval-upper-bound interval 3)))))

	    (else

	     (let ((multi-index (cons (i (cons j (cons k (cons l multi-index-tail)))))))

	       (or (and (not (every exact-integer? multi-index))

			(apply error "Interval-contains-multi-index?: at least one multi-index component is not an exact integer: " multi-index))

		   (and (not (= (##Interval-dimension interval) (length multi-index)))

			(apply error "Interval-contains-multi-index?: dimension of interval does not match number of arguments: " interval multi-index))

		   (every (lambda (l i u)

			    (and (<= l i)

				 (<  i u)))

			  (Interval-lower-bounds->list interval)

			  multi-index

			  (Interval-upper-bounds->list interval))))))))

;;; Applies f to every element of the domain, in lexicographical order

(define (Interval-for-each f interval)

  (or (and (not (Interval? interval))

	   (error "Interval-for-each: Argument is not a interval: " interval))

      (and (not (procedure? f))

	   (error "Interval-for-each: Argument is not a procedure: " f))

      (case (##Interval-dimension interval)

	((1) (let ((lower-i (##Interval-lower-bound interval 0))

		   (upper-i (##Interval-upper-bound interval 0)))

	       (let i-loop ((i lower-i))

		 (if (< i upper-i)

		     (begin

		       (f i)

		       (i-loop (+ i 1)))))))

	((2) (let ((lower-i (##Interval-lower-bound interval 0))

		   (lower-j (##Interval-lower-bound interval 1))

		   (upper-i (##Interval-upper-bound interval 0))

		   (upper-j (##Interval-upper-bound interval 1)))

	       (let i-loop ((i lower-i))

		 (if (< i upper-i)

		     (let j-loop ((j lower-j))

		       (if (< j upper-j)

			   (begin

			     (f i j)

			     (j-loop (+ j 1)))

			   (i-loop (+ i 1))))))))

	((3) (let ((lower-i (##Interval-lower-bound interval 0))

		   (lower-j (##Interval-lower-bound interval 1))

		   (lower-k (##Interval-lower-bound interval 2))

		   (upper-i (##Interval-upper-bound interval 0))

		   (upper-j (##Interval-upper-bound interval 1))

		   (upper-k (##Interval-upper-bound interval 2)))

	       (let i-loop ((i lower-i))

		 (if (< i upper-i)

		     (let j-loop ((j lower-j))

		       (if (< j upper-j)

			   (let k-loop ((k lower-k))

			     (if (< k upper-k)

				 (begin

				   (f i j k)

				   (k-loop (+ k 1)))

				 (j-loop (+ j 1))))

			   (i-loop (+ i 1))))))))

	((4) (let ((lower-i (##Interval-lower-bound interval 0))

		   (lower-j (##Interval-lower-bound interval 1))

		   (lower-k (##Interval-lower-bound interval 2))

		   (lower-l (##Interval-lower-bound interval 3))

		   (upper-i (##Interval-upper-bound interval 0))

		   (upper-j (##Interval-upper-bound interval 1))

		   (upper-k (##Interval-upper-bound interval 2))

		   (upper-l (##Interval-upper-bound interval 3)))

	       (let i-loop ((i lower-i))

		 (if (< i upper-i)

		     (let j-loop ((j lower-j))

		       (if (< j upper-j)

			   (let k-loop ((k lower-k))

			     (if (< k upper-k)

				 (let l-loop ((l lower-l))

				   (if (< l upper-l)

				       (begin

					 (f i j k l)

					 (l-loop (+ l 1)))

				       (k-loop (+ k 1))))

				 (j-loop (+ j 1))))

			   (i-loop (+ i 1))))))))

	(else

	 (let ((lower-bounds (list->vector (##Interval-lower-bounds->list interval)))

	       (upper-bounds (list->vector (##Interval-upper-bounds->list interval))))

	   (define (multi-index->args multi-index)

	     (cdr multi-index))

	   (define (vector->multi-index arg)

	     (define (tail-ref l i)

	       (let loop ((l l) (i i))

		 (if (= i 0)

		     l

		     (loop (cdr l) (- i 1)))))

	     (let ((l (vector->list arg))

		   (v (make-vector (vector-length arg))))

	       (do ((i (- (vector-length v) 1) (- i 1)))

		   ((< i 0) (cons v l))

		 (vector-set! v i (tail-ref l i)))))

	   (define (next-multi-index! multi-index lower-bounds upper-bounds)

	     (let ((v (car multi-index)))

	       (let loop ((i (- (vector-length v) 1)))

		 (and (<= 0 i)

		      (let* ((l-tail (vector-ref v i))

			     (next-index (+ (car l-tail) 1)))

			(if (< next-index (vector-ref upper-bounds i))

			    (begin

			      (set-car! l-tail next-index)

			      multi-index)

			    (begin

			      (set-car! l-tail (vector-ref lower-bounds i))

			      (loop (- i 1)))))))))

	   (let loop ((multi-index (vector->multi-index lower-bounds)))

	     (if multi-index

		 (begin

		   (apply f (multi-index->args multi-index))

		   (loop (next-multi-index! multi-index lower-bounds upper-bounds))))))))))

;;; I defined this separately in the hope that it would be speeded up.

;;; Calculates (...(operator (operator (operator identity (f multi-index_1)) (f multi-index_2)) (f multi-index_3)) ...)

;;; again in lexicographical order

(define (Interval-reduce f operator identity interval)

  (or (and (not (Interval? interval))

	   (error "Interval-reduce: Argument is not a interval: " interval))

      (and (not (procedure? f))

	   (error "Interval-reduce: Argument is not a procedure: " f))

      (and (not (procedure? operator))

	   (error "Interval-reduce: Operator is not a procedure: " operator))

      (case (##Interval-dimension interval)

	((1) (let ((lower-i (##Interval-lower-bound interval 0))

		   (upper-i (##Interval-upper-bound interval 0)))

	       (let i-loop ((i lower-i) (result identity))

		 (if (= i upper-i)

		     result

		     (i-loop (+ i 1) (operator result (f i)))))))

	((2) (let ((lower-i (##Interval-lower-bound interval 0))

		   (lower-j (##Interval-lower-bound interval 1))

		   (upper-i (##Interval-upper-bound interval 0))

		   (upper-j (##Interval-upper-bound interval 1)))

	       (let i-loop ((i lower-i) (result identity))

		 (if (= i upper-i)

		     result

		     (let j-loop ((j lower-j) (result result))

		       (if (= j upper-j)

			   (i-loop (+ i 1) result)

			   (j-loop (+ j 1) (operator result (f i j)))))))))

	((3) (let ((lower-i (##Interval-lower-bound interval 0))

		   (lower-j (##Interval-lower-bound interval 1))

		   (lower-k (##Interval-lower-bound interval 2))

		   (upper-i (##Interval-upper-bound interval 0))

		   (upper-j (##Interval-upper-bound interval 1))

		   (upper-k (##Interval-upper-bound interval 2)))

	       (let i-loop ((i lower-i) (result identity))

		 (if (= i upper-i)

		     result

		     (let j-loop ((j lower-j) (result result))

		       (if (= j upper-j)

			   (i-loop (+ i 1) result)

			   (let k-loop ((k lower-k) (result result))

			     (if (= k upper-k)

				 (j-loop (+ j 1) result)

				 (k-loop (+ k 1) (operator result (f i j k)))))))))))

	((4) (let ((lower-i (##Interval-lower-bound interval 0))

		   (lower-j (##Interval-lower-bound interval 1))

		   (lower-k (##Interval-lower-bound interval 2))

		   (lower-l (##Interval-lower-bound interval 3))

		   (upper-i (##Interval-upper-bound interval 0))

		   (upper-j (##Interval-upper-bound interval 1))

		   (upper-k (##Interval-upper-bound interval 2))

		   (upper-l (##Interval-upper-bound interval 3)))

	       (let i-loop ((i lower-i) (result identity))

		 (if (= i upper-i)

		     result

		     (let j-loop ((j lower-j) (result result))

		       (if (= j upper-j)

			   (i-loop (+ i 1) result)

			   (let k-loop ((k lower-k) (result result))

			     (if (= k upper-k)

				 (j-loop (+ j 1) result)

				 (let l-loop ((l lower-l) (result result))

				   (if (= l upper-l)

				       (k-loop (+ k 1) result)

				       (l-loop (+ l 1) (operator result (f i j k l)))))))))))))

	(else

	 (let ((lower-bounds (list->vector (##Interval-lower-bounds->list interval)))

	       (upper-bounds (list->vector (##Interval-upper-bounds->list interval))))

	   (define (tail-ref l i)

	     (let loop ((l l) (i i))

	       (if (= i 0)

		   l

		   (loop (cdr l) (- i 1)))))

	   (define (vector->multi-index arg)

	     (let ((l (vector->list arg))

		   (v (make-vector (vector-length arg))))

	       (do ((i (- (vector-length v) 1) (- i 1)))

		   ((< i 0) (cons v l))

		 (vector-set! v i (tail-ref l i)))))

	   (define (multi-index->args multi-index)

	     (cdr multi-index))

	   (define (next-multi-index! multi-index lower-bounds upper-bounds)

	     (let ((v (car multi-index)))

	       (let loop ((i (- (vector-length v) 1)))

		 (and (<= 0 i)

		      (let* ((l-tail (vector-ref v i))

			     (next-index (+ (car l-tail) 1)))

			(if (< next-index (vector-ref upper-bounds i))

			    (begin

			      (set-car! l-tail next-index)

			      multi-index)

			    (begin

			      (set-car! l-tail (vector-ref lower-bounds i))

			      (loop (- i 1)))))))))

	   (let loop ((multi-index (vector->multi-index lower-bounds))

		      (result identity))

	     (if multi-index

		 (let ((next-result (operator result (apply f (multi-index->args multi-index)))))

		   (loop (next-multi-index! multi-index lower-bounds upper-bounds)

			 next-result))

		 result)))))))

;; We'll use the same basic container for all types of arrays.

(define-structure Array-base

  ;; Part of all arrays

  domain                  ;; an Interval

  accessor                ;; (lambda (i_0 ... i_n-1) ...) returns a value for (i_0,...,i_n-1) in (Array-domain a)

  ;; Part of mutable arrays

  setter                  ;; (lambda (v i_0 ... i_n-1) ...) sets a value for (i_0,...,i_n-1) in (Array-domain a)

  ;; Part of Fixed arrays

  manipulators            ;; a Fixed-array-manipulator

  body                    ;; the backing store for this array

  indexer                 ;; see below

  affine?                 ;; whether the indexer is an affine map

  )

;; An array has a domain (which is an interval) and an accessor that maps that domain into some type of

;; Scheme objects

(define (build-Array domain accessor)

  (or (and (not (Interval? domain))

	   (error "build-Array: domain is not an interval: " domain))

      (and (not (procedure? accessor))

	   (error "build-Array: accessor is not a procedure: " accessor))

      (make-Array-base domain

		       accessor

		       #f        ; setter

		       #f        ; manipulators

		       #f        ; body

		       #f        ; indexer

		       #f        ; affine?

		       )))

(define (Array? obj)

  (Array-base? obj))

(define (Array-domain obj)

  (or (and (not (Array? obj))

	   (error "Array-domain: object is not an array: " obj))

      (Array-base-domain obj)))

(define (Array-accessor obj)

  (or (and (not (Array? obj))

	   (error "Array-accessor: object is not an array: " obj))

      (Array-base-accessor obj)))

(define (Array-lazy-curry array left-dimension)

  (or (and (not (Array? array))

	   (error "Array-lazy-curry: argument is not an array: " array left-dimension))

      (and (not (exact-integer? left-dimension))

	   (error "Array-lazy-curry: argument is not an exact integer: " array left-dimension))

      (and (not (< 0 left-dimension (##Interval-dimension (Array-domain array))))

	   (error "Array-lazy-curry: argument is not between 0 and (Interval-dimension (Array-domain array)) (exclusive): " array left-dimension))

      (call-with-values

	  (lambda () (Interval-curry (Array-domain array) left-dimension))

	(lambda (left-interval right-interval)

	  (let ((accessor (Array-accessor array)))

	    (case (##Interval-dimension left-interval)

	      ((1)  (case (##Interval-dimension right-interval)

		      ((1)  (build-Array left-interval

					 (lambda (i)      (build-Array right-interval (lambda (j)     (accessor i j))))))

		      ((2)  (build-Array left-interval

					 (lambda (i)      (build-Array right-interval (lambda (j k)   (accessor i j k))))))

		      ((3)  (build-Array left-interval

					 (lambda (i)      (build-Array right-interval (lambda (j k l) (accessor i j k l))))))

		      (else (build-Array left-interval

					 (lambda (i)      (build-Array right-interval (lambda multi-index (apply accessor i multi-index))))))))

	      ((2)  (case (##Interval-dimension right-interval)

		      ((1)  (build-Array left-interval

					 (lambda (i j)    (build-Array right-interval (lambda   (k)   (accessor i j k))))))

		      ((2)  (build-Array left-interval

					 (lambda (i j)    (build-Array right-interval (lambda   (k l) (accessor i j k l))))))

		      (else (build-Array left-interval

					 (lambda (i j)    (build-Array right-interval (lambda multi-index (apply accessor i j multi-index))))))))

	      ((3)  (case (##Interval-dimension right-interval)

		      ((1)  (build-Array left-interval

					 (lambda (i j k)  (build-Array right-interval (lambda     (l) (accessor i j k l))))))

		      (else (build-Array left-interval

					 (lambda (i j k)  (build-Array right-interval (lambda multi-index (apply accessor i j k multi-index))))))))

	      (else (build-Array left-interval

				 (lambda (#!rest left-multi-index)

				   (build-Array right-interval

						(lambda (#!rest right-multi-index)

						  (apply accessor (append left-multi-index right-multi-index)))))))))))))

(define (insert-arg-into-arg-list arg arg-list index)

  (define (helper arg-list i)

    (if (= i 0)

	(cons arg arg-list)

	(cons arg (helper (cdr arg-list) (- i 1)))))

  (helper arg-list index))

(define (Array-lazy-distinguish-one array index)

    (or (and (not (Array? array))

	   (error "Array-lazy-distinguish-one: argument is not an array: " array index))

      (and (not (exact-integer? index))

	   (error "Array-lazy-distinguish-one: argument is not an exact integer: " array index))

      (and (not (< 1 (##Interval-dimension (Array-domain array))))

	   (error "Array-lazy-distinguish-one:  (Interval-dimension (Array-domain array)) is not greater than 1 : " array index))

      (and (not (< -1 index (##Interval-dimension (Array-domain array))))

	   (error "Array-lazy-distinguish-one: argument is not between 0 (inclusive) and (Interval-dimension (Array-domain array)) (exclusive): " array index))

      (call-with-values

	  (lambda () (Interval-distinguish-one (Array-domain array) index))

	(lambda (outer-interval inner-interval)

	  (let ((accessor (Array-accessor array)))

	    (case (##Interval-dimension outer-interval)

	      ((1) (case index

		     ((0) (build-Array outer-interval

				       (lambda (j)     (build-Array inner-interval (lambda (i) (accessor i j))))))

		     ((1) (build-Array outer-interval

				       (lambda (i)     (build-Array inner-interval (lambda (j) (accessor i j))))))))

	      ((2) (case index

		     ((0) (build-Array outer-interval

				       (lambda (j k)   (build-Array inner-interval (lambda (i) (accessor i j k))))))

		     ((1) (build-Array outer-interval

				       (lambda (i k)   (build-Array inner-interval (lambda (j) (accessor i j k))))))

		     ((2) (build-Array outer-interval

				       (lambda (i j)   (build-Array inner-interval (lambda (k) (accessor i j k))))))))

	      ((3) (case index

		     ((0) (build-Array outer-interval

				       (lambda (j k l) (build-Array inner-interval (lambda (i) (accessor i j k l))))))

		     ((1) (build-Array outer-interval

				       (lambda (i k l) (build-Array inner-interval (lambda (j) (accessor i j k l))))))

		     ((2) (build-Array outer-interval

				       (lambda (i j l) (build-Array inner-interval (lambda (k) (accessor i j k l))))))

		     ((3) (build-Array outer-interval

				       (lambda (i j k) (build-Array inner-interval (lambda (l) (accessor i j k l))))))))

	      (else (build-Array outer-interval

				 (lambda (#!rest outer-index)

				   (build-Array inner-interval

						(lambda (m)

						  (apply accessor (insert-arg-into-arg-list m outer-index index)))))))))))))

#|

A Mutable array has, in addition a setter, that satisfies, roughly

If (i_1, ..., i_n)\neq (j_1, ..., j_n) \in (Array-domain a)

and

((Array-accessor a) j_1 ... j_n) => x

then "after"

((Array-setter a) v i_1 ... i_n)

we have

((Array-accessor a) j_1 ... j_n) => x

and

((Array-accessor a) i_1 ... i_n) => v

|#

(define (build-Mutable-array domain accessor setter)

  (or (and (not (Interval? domain))

	   (error "build-Mutable-array: domain is not an interval: " domain))

      (and (not (procedure? accessor))

	   (error "build-Mutable-array: accessor is not a procedure: " accessor))

      (and (not (procedure? setter))

	   (error "build-Mutable-array: setter is not a procedure: " setter))

      (make-Array-base domain

		       accessor

		       setter

		       #f        ; manipulators

		       #f        ; body

		       #f        ; indexer

		       #f        ; affine?

		       )))

(define (Mutable-array? obj)

  (and (Array? obj)

       (not (eq? (Array-base-setter obj) #f))))

(define (Mutable-array-domain obj)

  (or (and (not (Mutable-array? obj))

	   (error "Mutable-array-domain: object is not an mutable array: " obj))

      (Array-base-domain obj)))

(define (Mutable-array-accessor obj)

  (or (and (not (Mutable-array? obj))

	   (error "Mutable-array-accessor: object is not an mutable array: " obj))

      (Array-base-accessor obj)))

(define (Mutable-array-setter obj)

  (or (and (not (Mutable-array? obj))

	   (error "Mutable-array-setter: object is not an mutable array: " obj))

      (Array-base-setter obj)))

(define (Mutable-array-lazy-curry array left-dimension)

  (declare (fixnum))

  (or (and (not (Array? array))

	   (error "Mutable-array-lazy-curry: argument is not a Mutable-array: " array left-dimension))

      (and (not (exact-integer? left-dimension))

	   (error "Mutable-array-lazy-curry: argument is not an exact integer: " array left-dimension))

      (and (not (< 0 left-dimension (##Interval-dimension (Mutable-array-domain array))))

	   (error "Mutable-array-lazy-curry: argument is not between 0 and (Interval-dimension (Mutable-array-domain array)) (exclusive): " array left-dimension))

      (call-with-values

	  (lambda () (Interval-curry (Mutable-array-domain array) left-dimension))

	(lambda (left-interval right-interval)

	  (let ((accessor (Mutable-array-accessor array))

		(setter   (Mutable-array-setter   array)))

	    (case (##Interval-dimension left-interval)

	      ((1)  (case (##Interval-dimension right-interval)

		      ((1)  (build-Array left-interval

					 (lambda (i)

					   (build-Mutable-array right-interval

							(lambda (j)

							  (accessor i j))

							(lambda (v j)

							  (setter v i j))))))

		      ((2)  (build-Array left-interval

					 (lambda (i)

					   (build-Mutable-array right-interval

							(lambda (j k)

							  (accessor i j k))

							(lambda (v j k)

							  (setter v i j k))))))

		      ((3)  (build-Array left-interval

					 (lambda (i)

					   (build-Mutable-array right-interval

							(lambda (j k l)

							  (accessor i j k l))

							(lambda (v j k l)

							  (setter v i j k l))))))

		      (else (build-Array left-interval

					 (lambda (i)

					   (build-Mutable-array right-interval

							(lambda (#!rest multi-index)

							  (apply accessor i multi-index))

							(lambda (v #!rest multi-index)

							  (apply setter v i multi-index))))))))

	      ((2)  (case (##Interval-dimension right-interval)

		      ((1)  (build-Array left-interval

					 (lambda (i j)

					   (build-Mutable-array right-interval

							(lambda (k)

							  (accessor i j k))

							(lambda (v k)

							  (setter v i j k))))))

		      ((2)  (build-Array left-interval

					 (lambda (i j)

					   (build-Mutable-array right-interval

							(lambda (k l)

							  (accessor i j k l))

							(lambda (v k l)

							  (setter v i j k l))))))

		      (else (build-Array left-interval

					 (lambda (i j)

					   (build-Mutable-array right-interval

							(lambda (#!rest multi-index)

							  (apply accessor i j multi-index))

							(lambda (v #!rest multi-index)

							  (apply setter v i j multi-index))))))))

	      ((3)  (case (##Interval-dimension right-interval)

		      ((1)  (build-Array left-interval

					 (lambda (i j k)

					   (build-Mutable-array right-interval

							(lambda (l)

							  (accessor i j k l))

							(lambda (v l)

							  (setter v i j k l))))))

		      (else (build-Array left-interval

					 (lambda (i j k)

					   (build-Mutable-array right-interval

							(lambda (#!rest multi-index)

							  (apply accessor i j k multi-index))

							(lambda (v #!rest multi-index)

							  (apply setter v i j k multi-index))))))))

	      (else (build-Array left-interval

				 (lambda (#!rest left-multi-index)

				   (build-Mutable-array right-interval

						(lambda (#!rest right-multi-index)

						  (apply accessor (append left-multi-index right-multi-index)))

						(lambda (v #!rest right-multi-index)

						  (apply setter v (append left-multi-index right-multi-index)))))))))))))

(define (Mutable-array-lazy-distinguish-one array index)

    (or (and (not (Mutable-array? array))

	   (error "Mutable-array-lazy-distinguish-one: argument is not an array: " array index))

      (and (not (exact-integer? index))

	   (error "Mutable-array-lazy-distinguish-one: argument is not an exact integer: " array index))

      (and (not (< 1 (##Interval-dimension (Array-domain array))))

	   (error "Mutable-array-lazy-distinguish-one:  (Interval-dimension (Array-domain array)) is not greater than 1 : " array index))

      (and (not (< -1 index (##Interval-dimension (Array-domain array))))

	   (error "Mutable-array-lazy-distinguish-one: argument is not between 0 (inclusive) and (Interval-dimension (Array-domain array)) (exclusive): " array index))

      (call-with-values

	  (lambda () (Interval-distinguish-one (Array-domain array) index))

	(lambda (outer-interval inner-interval)

	  (let ((accessor (Mutable-array-accessor array))

		(setter   (Mutable-array-setter array)))

	    (case (##Interval-dimension outer-interval)

	      ((1) (case index

		     ((0) (build-Array outer-interval

				       (lambda (j)     (build-Mutable-array inner-interval

									    (lambda (  i) (accessor i j))

									    (lambda (v i) (setter v i j))))))

		     ((1) (build-Array outer-interval

				       (lambda (i)     (build-Mutable-array inner-interval

									    (lambda (  j) (accessor i j))

									    (lambda (v j) (setter v i j))))))))

	      ((2) (case index

		     ((0) (build-Array outer-interval

				       (lambda (j k)   (build-Mutable-array inner-interval

									    (lambda (  i) (accessor i j k))

									    (lambda (v i) (setter v i j k))))))

		     ((1) (build-Array outer-interval

				       (lambda (i k)   (build-Mutable-array inner-interval

									    (lambda (  j) (accessor i j k))

									    (lambda (v j) (setter v i j k))))))

		     ((2) (build-Array outer-interval

				       (lambda (i j)   (build-Mutable-array inner-interval

									    (lambda (  k) (accessor i j k))

									    (lambda (v k) (setter v i j k))))))))

	      ((3) (case index

		     ((0) (build-Array outer-interval

				       (lambda (j k l) (build-Mutable-array inner-interval

									    (lambda (  i) (accessor i j k l))

									    (lambda (v i) (setter v i j k l))))))

		     ((1) (build-Array outer-interval

				       (lambda (i k l) (build-Mutable-array inner-interval

									    (lambda (  j) (accessor i j k l))

									    (lambda (v j) (setter v i j k l))))))

		     ((2) (build-Array outer-interval

				       (lambda (i j l) (build-Mutable-array inner-interval

									    (lambda (  k) (accessor i j k l))

									    (lambda (v k) (setter v i j k l))))))

		     ((3) (build-Array outer-interval

				       (lambda (i j k) (build-Mutable-array inner-interval

									    (lambda (  l) (accessor i j k l))

									    (lambda (v l) (setter v i j k l))))))))

	      (else (build-Array outer-interval

				 (lambda (#!rest outer-index)

				   (build-Array inner-interval

						(lambda (  m)

						  (apply accessor (insert-arg-into-arg-list m outer-index index)))

						(lambda (v m)

						  (apply setter v (insert-arg-into-arg-list m outer-index index)))))))))))))

#|

Fixed-array-manipulators contains functions and objects to manipulate the

backing store of a Fixed-array.

accessor: (lambda (body i) ...)   returns the value of body at index i

setter:   (lambda (body i v) ...) sets the value of body at index i to v

maker:    (lambda (n) ...)        makes a body of length n

length:   (lambda (body) ...)     returns the number of objects in body

default:  object                  is the default value with which to fill body

|#

(define-structure Fixed-array-manipulators accessor setter maker length default)

#|

We define specialized array manipulators for:

32- and 64-bit floating-point numbers,

8-, 16-, 32-, and 64-bit signed integers,

8-, 16-, 32-, and 64-bit unsigned integers, and

1-bit unsigned integers

as well as generic objects.

|#

(define-macro (make-array-manipulators)

  (define (symbol-concatenate . symbols)

    (string->symbol (apply string-append (map (lambda (s)

						(if (string? s)

						    s

						    (symbol->string s)))

					      symbols))))

  `(begin

     ,@(map (lambda (name prefix default)

	      `(define ,(symbol-concatenate name '-array-manipulators)

		 (make-Fixed-array-manipulators

		  ;; accessor:

		  (lambda (v i)

		    (,(symbol-concatenate prefix 'vector-ref) v i))

		  ;; setter:

		  (lambda (v i val)

		    (,(symbol-concatenate prefix 'vector-set!) v i val))

		  ;; maker:

		  ,(symbol-concatenate 'make- prefix 'vector)

		  ;; length:

		  ,(symbol-concatenate prefix 'vector-length)

		  ;; default:

		  ,default)))

	    '(generic s8 u8 s16 u16 s32 u32 s64 u64 f32 f64)

	    '(""      s8 u8 s16 u16 s32 u32 s64 u64 f32 f64)

	    '(#f       0  0   0   0   0   0   0   0 0.0 0.0))))

(make-array-manipulators)

;;; for bit-arrays, body is a vector, the first element of which is the actual number of elements,

;;; the second element of which is a u16vector that contains the bit string

(define u1-array-manipulators

  (make-Fixed-array-manipulators

   ;; accessor:

   (lambda (v i)

     (let ((index (##fixnum.arithmetic-shift-right i 4))

	   (shift (##fixnum.bitwise-and i 15))

	   (bodyv (vector-ref v  1)))

       (##fixnum.bitwise-and

	(##fixnum.arithmetic-shift-right

	 (u16vector-ref bodyv index)

	 shift)

	1)))

   ;; setter: 

   (lambda (v i val)

     (let ((index (##fixnum.arithmetic-shift-right i 4))

	   (shift (##fixnum.bitwise-and i 15))

	   (bodyv (vector-ref v  1)))

       (u16vector-set! bodyv index (##fixnum.bitwise-ior

				    (##fixnum.arithmetic-shift-left val shift)

				    (##fixnum.bitwise-and

				     (u16vector-ref bodyv index)

				     (##fixnum.bitwise-not

				      (##fixnum.arithmetic-shift-left 1 shift)))))))

   ;; maker:

   (lambda (size initializer)

     (let ((u16-size (##fixnum.arithmetic-shift-right (+ size 15) 4)))

       (vector size (make-u16vector u16-size (if (zero? initializer) 0 65535)))))

   ;; length:

   (lambda (v)

     (vector-ref v 0))

   ;; default:

   0))

#|

Conceptually, an indexer is itself a 1-1 Array that maps one interval to another; thus, it is

an example of an array that can return multiple values.

Rather than trying to formalize this idea, and trying to get it to work with Array-map,

Array-reduce, ..., we'll just manipulate the accessor functions of these conceptual Arrays.

Affine indexers are 1-1 affine maps from one interval to another.

The indexer field of a Fixed-array obj is a 1-1 mapping from

(Array-domain obj)

to [0, top), where top is 

((Fixed-array-manipulators-length (Fixed-array-manipulators obj)) (Fixed-array-body obj))

|#

#|

build-affine-indexer specializes

(define (my-build-indexer base lower-bounds increments)

  (lambda multi-index

    (apply + base (map * increments (map - multi-index lower-bounds)))))

|#

(define (build-affine-indexer base lower-bounds increments)

  (case (length lower-bounds)

    ((1)  (build-indexer-1 base

			   (car lower-bounds)

			   (car increments)))

    ((2)  (build-indexer-2 base

			   (car lower-bounds) (cadr lower-bounds)

			   (car increments)   (cadr increments)))

    ((3)  (build-indexer-3 base

			   (car lower-bounds) (cadr lower-bounds) (caddr lower-bounds)

			   (car increments)   (cadr increments)   (caddr increments)))

    ((4)  (build-indexer-4 base

			   (car lower-bounds) (cadr lower-bounds) (caddr lower-bounds) (cadddr lower-bounds)

			   (car increments)   (cadr increments)   (caddr increments)   (cadddr increments)))

    (else (build-indexer-generic base lower-bounds increments))))

;; unfortunately, the next two functions were written by hand, so beware of bugs.

(define (build-indexer-1 base

			 low-0

			 increment-0)

  (if (zero? base)

      (if (zero? low-0)

	  (cond ((= 1 increment-0)    (lambda (i) i))

		((= -1 increment-0)   (lambda (i) (- i)))

		(else                 (lambda (i) (* i increment-0))))

	  (cond ((= 1 increment-0)    (lambda (i) (- i low-0)))

		((= -1 increment-0)   (lambda (i) (- low-0 i)))

		(else                 (lambda (i) (* increment-0 (- i low-0))))))

      (if (zero? low-0)

	  (cond ((= 1 increment-0)    (lambda (i) (+ base i)))

		((= -1 increment-0)   (lambda (i) (- base i)))

		(else                 (lambda (i) (+ base (* increment-0 i)))))

	  (cond ((= 1 increment-0)    (lambda (i) (+ base (- i low-0))))

		((= -1 increment-0)   (lambda (i) (+ base (- low-0 i))))

		(else                 (lambda (i) (+ base (* increment-0 (- i low-0)))))))))

(define (build-indexer-2 base

			 low-0       low-1

			 increment-0 increment-1)

  (if (zero? base)

      (if (zero? low-0)

	  (cond ((= 1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)    (lambda (i j) (+ i j)))

			   ((= -1 increment-1)   (lambda (i j) (+ i (- j))))

			   (else                 (lambda (i j) (+ i (* increment-1 j)))))

		     (cond ((= 1 increment-1)    (lambda (i j) (+ i (- j low-1))))

			   ((= -1 increment-1)   (lambda (i j) (+ i (- low-1 j))))

			   (else                 (lambda (i j) (+ i (* increment-1 (- j low-1))))))))

		((= -1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (- j                           i)))

			   ((= -1 increment-1)  (lambda (i j) (- (- j)                       i)))

			   (else                (lambda (i j) (- (* increment-1 j)           i))))

		     (cond ((= 1 increment-1)   (lambda (i j) (- (- j low-1)                 i)))

			   ((= -1 increment-1)  (lambda (i j) (- (- low-1 j)                 i)))

			   (else                (lambda (i j) (- (* increment-1 (- j low-1)) i))))))

		(else

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (+ (* increment-0 i) j)))

			   ((= -1 increment-1)  (lambda (i j) (+ (* increment-0 i) (- j))))

			   (else                (lambda (i j) (+ (* increment-0 i) (* increment-1 j)))))

		     (cond ((= 1 increment-1)   (lambda (i j) (+ (* increment-0 i) (- j low-1))))

			   ((= -1 increment-1)  (lambda (i j) (+ (* increment-0 i) (- low-1 j))))

			   (else                (lambda (i j) (+ (* increment-0 i) (* increment-1 (- j low-1)))))))))

	  (cond ((= 1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)    (lambda (i j) (+ (- i low-0) j)))

			   ((= -1 increment-1)   (lambda (i j) (+ (- i low-0) (- j))))

			   (else                 (lambda (i j) (+ (- i low-0) (* increment-1 j)))))

		     (cond ((= 1 increment-1)    (lambda (i j) (+ (- i low-0) (- j low-1))))

			   ((= -1 increment-1)   (lambda (i j) (+ (- i low-0) (- low-1 j))))

			   (else                 (lambda (i j) (+ (- i low-0) (* increment-1 (- j low-1))))))))

		((= -1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (- j                           (- i low-0))))

			   ((= -1 increment-1)  (lambda (i j) (- (- j)                       (- i low-0))))

			   (else                (lambda (i j) (- (* increment-1 j)           (- i low-0)))))

		     (cond ((= 1 increment-1)   (lambda (i j) (- (- j low-1)                 (- i low-0))))

			   ((= -1 increment-1)  (lambda (i j) (- (- low-1 j)                 (- i low-0))))

			   (else                (lambda (i j) (- (* increment-1 (- j low-1)) (- i low-0)))))))

		(else

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (+ (* increment-0 (- i low-0)) j)))

			   ((= -1 increment-1)  (lambda (i j) (+ (* increment-0 (- i low-0)) (- j))))

			   (else                (lambda (i j) (+ (* increment-0 (- i low-0)) (* increment-1 j)))))

		     (cond ((= 1 increment-1)   (lambda (i j) (+ (* increment-0 (- i low-0)) (- j low-1))))

			   ((= -1 increment-1)  (lambda (i j) (+ (* increment-0 (- i low-0)) (- low-1 j))))

			   (else                (lambda (i j) (+ (* increment-0 (- i low-0)) (* increment-1 (- j low-1))))))))))

      (if (zero? low-0)

	  (cond ((= 1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)    (lambda (i j) (+ base i j)))

			   ((= -1 increment-1)   (lambda (i j) (+ base i (- j))))

			   (else                 (lambda (i j) (+ base i (* increment-1 j)))))

		     (cond ((= 1 increment-1)    (lambda (i j) (+ base i (- j low-1))))

			   ((= -1 increment-1)   (lambda (i j) (+ base i (- low-1 j))))

			   (else                 (lambda (i j) (+ base i (* increment-1 (- j low-1))))))))

		((= -1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (- (+ base j)                           i)))

			   ((= -1 increment-1)  (lambda (i j) (- (- base j)                           i)))

			   (else                (lambda (i j) (- (+ base (* increment-1 j))           i))))

		     (cond ((= 1 increment-1)   (lambda (i j) (- (+ base (- j low-1))                 i)))

			   ((= -1 increment-1)  (lambda (i j) (- (+ base (- low-1 j))                 i)))

			   (else                (lambda (i j) (- (+ base (* increment-1 (- j low-1))) i))))))

		(else

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (+ base (* increment-0 i) j)))

			   ((= -1 increment-1)  (lambda (i j) (+ base (* increment-0 i) (- j))))

			   (else                (lambda (i j) (+ base (* increment-0 i) (* increment-1 j)))))

		     (cond ((= 1 increment-1)   (lambda (i j) (+ base (* increment-0 i) (- j low-1))))

			   ((= -1 increment-1)  (lambda (i j) (+ base (* increment-0 i) (- low-1 j))))

			   (else                (lambda (i j) (+ base (* increment-0 i) (* increment-1 (- j low-1)))))))))

	  (cond ((= 1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)    (lambda (i j) (+ base (- i low-0) j)))

			   ((= -1 increment-1)   (lambda (i j) (+ base (- i low-0) (- j))))

			   (else                 (lambda (i j) (+ base (- i low-0) (* increment-1 j)))))

		     (cond ((= 1 increment-1)    (lambda (i j) (+ base (- i low-0) (- j low-1))))

			   ((= -1 increment-1)   (lambda (i j) (+ base (- i low-0) (- low-1 j))))

			   (else                 (lambda (i j) (+ base (- i low-0) (* increment-1 (- j low-1))))))))

		((= -1 increment-0)

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (- (+ base j)                           (- i low-0))))

			   ((= -1 increment-1)  (lambda (i j) (- (- base j)                           (- i low-0))))

			   (else                (lambda (i j) (- (+ base (* increment-1 j))           (- i low-0)))))

		     (cond ((= 1 increment-1)   (lambda (i j) (- (+ base (- j low-1))                 (- i low-0))))

			   ((= -1 increment-1)  (lambda (i j) (- (+ base (- low-1 j))                 (- i low-0))))

			   (else                (lambda (i j) (- (+ base (* increment-1 (- j low-1))) (- i low-0)))))))

		(else

		 (if (zero? low-1)

		     (cond ((= 1 increment-1)   (lambda (i j) (+ base (* increment-0 (- i low-0)) j)))

			   ((= -1 increment-1)  (lambda (i j) (+ base (* increment-0 (- i low-0)) (- j))))

			   (else                (lambda (i j) (+ base (* increment-0 (- i low-0)) (* increment-1 j)))))

		     (cond ((= 1 increment-1)   (lambda (i j) (+ base (* increment-0 (- i low-0)) (- j low-1))))

			   ((= -1 increment-1)  (lambda (i j) (+ base (* increment-0 (- i low-0)) (- low-1 j))))

			   (else                (lambda (i j) (+ base (* increment-0 (- i low-0)) (* increment-1 (- j low-1))))))))))))

;;; after this we basically punt

(define (build-indexer-3 base

			 low-0       low-1       low-2

			 increment-0 increment-1 increment-2)

  (if (= 0 low-0 low-1 low-2)

      (if (= base 0)

	  (if (= increment-2 1)

	      (lambda (i j k)

		(+ (* increment-0 i)

		   (* increment-1 j)

		   k))

	      (lambda (i j k)

		(+ (* increment-0 i)

		   (* increment-1 j)

		   (* increment-2 k))))

	  (if (= increment-2 1)

	      (lambda (i j k)

		(+ base

		   (* increment-0 i)

		   (* increment-1 j)

		   k))

	      (lambda (i j k)

		(+ base

		   (* increment-0 i)

		   (* increment-1 j)

		   (* increment-2 k)))))

      (if (= base 0)

	  (if (= increment-2 1)

	      (lambda (i j k)

		(+ (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (- k low-2)))

	      (lambda (i j k)

		(+ (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (* increment-2 (- k low-2)))))

	  (if (= increment-2 1)

	      (lambda (i j k)

		(+ base

		   (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (- k low-2)))

	      (lambda (i j k)

		(+ base

		   (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (* increment-2 (- k low-2))))))))

(define (build-indexer-4 base

			 low-0       low-1       low-2       low-3

			 increment-0 increment-1 increment-2 increment-3)

  (if (= 0 low-0 low-1 low-2 low-3)

      (if (= base 0)

	  (if (= increment-3 1)

	      (lambda (i j k l)

		(+ (* increment-0 i)

		   (* increment-1 j)

		   (* increment-2 k)

		   l))

	      (lambda (i j k l)

		(+ (* increment-0 i)

		   (* increment-1 j)

		   (* increment-2 k)

		   (* increment-3 l))))

	  (if (= increment-3 1)

	      (lambda (i j k l)

		(+ base

		   (* increment-0 i)

		   (* increment-1 j)

		   (* increment-2 k)

		   l))

	      (lambda (i j k l)

		(+ base

		   (* increment-0 i)

		   (* increment-1 j)

		   (* increment-2 k)

		   (* increment-3 l)))))

      (if (= base 0)

	  (if (= increment-3 1)

	      (lambda (i j k l)

		(+ (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (* increment-2 (- k low-2))

		   (- l low-3)))

	      (lambda (i j k l)

		(+ (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (* increment-2 (- k low-2))

		   (* increment-3 (- l low-3)))))

	  (if (= increment-3 1)

	      (lambda (i j k l)

		(+ base

		   (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (* increment-2 (- k low-2))

		   (- l low-3)))

	      (lambda (i j k l)

		(+ base

		   (* increment-0 (- i low-0))

		   (* increment-1 (- j low-1))

		   (* increment-2 (- k low-2))

		   (* increment-3 (- l low-3))))))))

(define (build-indexer-generic base lower-bounds increments)

  (lambda multi-index

    (do ((multi-index  multi-index  (cdr multi-index))

	 (lower-bounds lower-bounds (cdr lower-bounds))

	 (increments   increments   (cdr increments))

	 (result       base         (+ result (* (car increments)

						 (- (car multi-index)

						    (car lower-bounds))))))

	((null? multi-index) result))))

(define (affine-indexer= indexer1 indexer2 interval)

  (or (and (not (Interval? interval))

	   (error "affine-indexer=: argument is not an interval " interval))

      (and (not (procedure? indexer1))

	   (error "affine-indexer=: argument is not a procedure " indexer1))

      (and (not (procedure? indexer2))

	   (error "affine-indexer=: argument is not a procedure " indexer2))

      (let ((lower-bounds (Interval-lower-bounds->list interval))

	    (upper-bounds (Interval-upper-bounds->list interval)))

	(let ((multi-index (map (lambda (x) x) lower-bounds)))

	  (let outer ((i multi-index)

		      (l lower-bounds)

		      (u upper-bounds))

	    (and (= (apply indexer1 multi-index)

		    (apply indexer2 multi-index))

		 (or (null? i)

		     (begin

		       (if (< (+ 1 (car l)) (car u))

			   (set-car! i (+ (car i) 1)))

		       (outer (cdr i) (cdr l) (cdr u))))))))))

(define (generic-indexer= indexer1 indexer2 interval)

  (or (and (not (Interval? interval))

	   (error "generic-indexer=: argument is not an interval " interval))

      (and (not (procedure? indexer1))

	   (error "generic-indexer=: argument is not a procedure " indexer1))

      (and (not (procedure? indexer2))

	   (error "generic-indexer=: argument is not a procedure " indexer2))

      (Interval-reduce (lambda multi-index

			 (= (apply indexer1 multi-index)

			    (apply indexer2 multi-index)))

		       (lambda (x y) (and x y))

		       #t

		       interval)))

#|

The default accessor and the setter of a Fixed-array a are given by

(lambda (i_0 ... i_n-1)

  ((Fixed-array-manipulators-accessor a)

   (Fixed-array-body a)

   ((Fixed-array-indexer a) i_0 ... i_n-1)))

(lambda (v i_0 ... i_n-1)

  ((Fixed-array-manipulators-setter a)

   (Fixed-array-body a)

   ((Fixed-array-indexer a) i_0 ... i_n-1)

   v))

The default initializer-value is 

(Fixed-array-manipulators-default manipulators)

The default body is

((Fixed-array-manipulators-maker manipulators)

 (Interval-volume domain)

 initializer-value)

The default indexer is the mapping of

the domain to the natural numbers in lexicographical order.

|#

(define (Fixed-array? obj)

  (and (Mutable-array? obj)

       (not (eq? (Array-base-body obj) #f))))

(define (Fixed-array-domain obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-domain: object is not an mutable array: " obj))

      (Array-base-domain obj)))

(define (Fixed-array-accessor obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-accessor: object is not an mutable array: " obj))

      (Array-base-accessor obj)))

(define (Fixed-array-setter obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-setter: object is not an mutable array: " obj))

      (Array-base-setter obj)))

(define (Fixed-array-body obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-body: argument is not a fixed array: " obj))

      (Array-base-body obj)))

(define (Fixed-array-indexer obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-indexer: argument is not a fixed array: " obj))

      (Array-base-indexer obj)))

(define (Fixed-array-manipulators obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-manipulators: argument is not a fixed array: " obj))

      (Array-base-manipulators obj)))

(define (Fixed-array-affine? obj)

  (or (and (not (Fixed-array? obj))

	   (error "Fixed-array-affine?: argument is not a fixed array: " obj))

      (Array-base-affine? obj)))

(define (##finish-Fixed-array domain manipulators body indexer affine?)

  (let ((accessor (let ((manipulator-accessor (Fixed-array-manipulators-accessor manipulators))

			(indexer indexer)

			(body body))

		    (case (##Interval-dimension domain)

		      ((1)  (lambda (i)         (manipulator-accessor body (indexer i))))

		      ((2)  (lambda (i j)       (manipulator-accessor body (indexer i j))))

		      ((3)  (lambda (i j k)     (manipulator-accessor body (indexer i j k))))

		      ((4)  (lambda (i j k l)   (manipulator-accessor body (indexer i j k l))))

		      (else (lambda multi-index (manipulator-accessor body (apply indexer multi-index)))))))

	(setter (let ((manipulator-setter (Fixed-array-manipulators-setter manipulators))

		      (indexer indexer)

		      (body body))

		  (case (##Interval-dimension domain)

		    ((1)  (lambda (value i)             (manipulator-setter body (indexer i)                 value)))

		    ((2)  (lambda (value i j)           (manipulator-setter body (indexer i j)               value)))

		    ((3)  (lambda (value i j k)         (manipulator-setter body (indexer i j k)             value)))

		    ((4)  (lambda (value i j k l)       (manipulator-setter body (indexer i j k l)           value)))

		    (else (lambda (value . multi-index) (manipulator-setter body (apply indexer multi-index) value)))))))

    (make-Array-base domain

		     accessor

		     setter

		     manipulators

		     body

		     indexer

		     affine?)))

(define (build-Fixed-array #!key

			   (domain            (macro-absent-obj))

			   (manipulators      (macro-absent-obj))

			   (body              (macro-absent-obj))

			   (indexer           (macro-absent-obj))

			   (initializer-value (macro-absent-obj))

			   (affine?           (macro-absent-obj)))

  (or (and (eq? (macro-absent-obj) domain)

	   (error "build-Fixed-array: the domain must be given: " domain))

      (and (not (Interval? domain))

	   (error "build-Fixed-array: domain is not an interval: " domain))

      (and (eq? (macro-absent-obj) manipulators)

	   (error "build-Fixed-array: the manipulators must be given: " manipulators))

      (and (not (Fixed-array-manipulators? manipulators))

	   (error "build-Fixed-array: manipulators are not fixed-array-manipulators: " manipulators))

      (and (not (eq? (macro-absent-obj) indexer))

	   (not (procedure? indexer))

	   (error "build-Fixed-array:  indexer is not a procedure: " indexer))

      (and (and (not (eq? affine? (macro-absent-obj)))

		(not (boolean? affine?)))

	   (error "build-Fixed-array: affine? must be a boolean: " affine?))

      (and (and (not (eq? body              (macro-absent-obj)))

		(not (eq? initializer-value (macro-absent-obj))))

	   (error "build-Fixed-array: you can not pass both body and initializer-value as arguments " body initializer-value))

      (let* ((affine? (if (boolean? affine?)

			  affine?

			  (eq? indexer (macro-absent-obj))))

	     (initializer-value (if (eq? initializer-value (macro-absent-obj))

				    (Fixed-array-manipulators-default manipulators)

				    initializer-value))

	     (body (if (eq? body (macro-absent-obj))

		       ((Fixed-array-manipulators-maker manipulators) (Interval-volume domain) initializer-value)

		       body))

	     (indexer (if (eq? indexer (macro-absent-obj))

			  (case (Interval-dimension domain)

			    ((1) (let ((low-0 (Interval-lower-bound domain 0))

				       (increment-0 1))

				   (build-indexer-1 0 low-0 increment-0)))

			    ((2) (let* ((low-0 (##Interval-lower-bound domain 0))

					(low-1 (##Interval-lower-bound domain 1))

					(increment-1 1)

					(increment-0 (* increment-1

							(- (##Interval-upper-bound domain 1)

							   (##Interval-lower-bound domain 1)))))

				   (build-indexer-2 0

						    low-0 low-1

						    increment-0 increment-1)))

			    ((3) (let* ((low-0 (##Interval-lower-bound domain 0))

					(low-1 (##Interval-lower-bound domain 1))

					(low-2 (##Interval-lower-bound domain 2))

					(increment-2 1)

					(increment-1 (* increment-2

							(- (##Interval-upper-bound domain 2)

							   (##Interval-lower-bound domain 2))))

					(increment-0 (* increment-1

							(- (##Interval-upper-bound domain 1)

							   (##Interval-lower-bound domain 1)))))

				   (build-indexer-3 0

						    low-0 low-1 low-2

						    increment-0 increment-1 increment-2)))

			    ((4) (let* ((low-0 (##Interval-lower-bound domain 0))

					(low-1 (##Interval-lower-bound domain 1))

					(low-2 (##Interval-lower-bound domain 2))

					(low-3 (##Interval-lower-bound domain 3))

					(increment-3 1)

					(increment-2 (* increment-3

							(- (##Interval-upper-bound domain 3)

							   (##Interval-lower-bound domain 3))))

					(increment-1 (* increment-2

							(- (##Interval-upper-bound domain 2)

							   (##Interval-lower-bound domain 2))))

					(increment-0 (* increment-1

							(- (##Interval-upper-bound domain 1)

							   (##Interval-lower-bound domain 1)))))

				   (build-indexer-4 0

						    low-0 low-1 low-2 low-3

						    increment-0 increment-1 increment-2 increment-3)))

			    (else

			     (let ((lower-bounds (##Interval-lower-bounds->list domain))

				   (upper-bounds (##Interval-upper-bounds->list domain)))

			       (let ((ranges (map (lambda (u l) (- u l)) upper-bounds lower-bounds)))

				 (do ((ranges (reverse ranges) (cdr ranges))

				      (increments (list 1) (cons (* (car increments) (car ranges))

								 increments)))

				     ((null? (cdr ranges)) (build-indexer-generic 0 lower-bounds increments)))))))

			  indexer)))

	(##finish-Fixed-array domain

			      manipulators

			      body

			      indexer

			      affine?))))

#|

The domain of the result is the same as the domain of the argument.

Builds a new Fixed-array and populates the body of the result in

lexicographical order of the domain.

|#

(define (Array->Fixed-array result-manipulators array)

  (or (and (not (Array? array))

	   (error "Array->Fixed-array: Argument is not an array: " array))

      (and (not (Fixed-array-manipulators? result-manipulators))

	   (error "Array->Fixed-array: Argument is not a set of fixed-array-manipulators: " result-manipulators))

      (let ((domain (Array-domain array)))

	(let ((result (build-Fixed-array

		       domain: domain

		       manipulators: result-manipulators)))

	  (let ((range-setter (Mutable-array-setter result))

		(domain-accessor (Array-accessor array)))

	    (Interval-for-each (case (##Interval-dimension domain)

				 ((1)  (lambda (i)       (range-setter (domain-accessor i)       i)))

				 ((2)  (lambda (i j)     (range-setter (domain-accessor i j)     i j)))

				 ((3)  (lambda (i j k)   (range-setter (domain-accessor i j k)   i j k)))

				 ((4)  (lambda (i j k l) (range-setter (domain-accessor i j k l) i j k l)))

				 (else (lambda multi-index

					 (apply range-setter (apply domain-accessor multi-index) multi-index))))

			       domain)

	    result)))))

#|

In the next function, old-indexer is an affine 1-1 mapping from an interval to [0,N), for some N.

new-domain->old-domain is an affine 1-1 mapping from new-domain to the domain of old-indexer.

|#

(define (compose-affine-indexers old-indexer new-domain new-domain->old-domain)

  (case (##Interval-dimension new-domain)

    ((1) (let* ((lower-0 (##Interval-lower-bound new-domain 0))

		(upper-0 (##Interval-upper-bound new-domain 0))

		(base (call-with-values

			  (lambda () (new-domain->old-domain lower-0))

			old-indexer))

		(increment-0 (if (< (+ lower-0 1) upper-0)

				 (- (call-with-values

					(lambda () (new-domain->old-domain (+ lower-0 1)))

				      old-indexer)

				    base)

				 0)))

	   (build-indexer-1 base lower-0 increment-0)))

    ((2) (let* ((lower-0 (##Interval-lower-bound new-domain 0))

		(lower-1 (##Interval-lower-bound new-domain 1))

		(upper-0 (##Interval-upper-bound new-domain 0))

		(upper-1 (##Interval-upper-bound new-domain 1))

		(base (call-with-values

			  (lambda () (new-domain->old-domain lower-0 lower-1))

			old-indexer))

		(increment-0 (if (< (+ lower-0 1) upper-0)

				 (- (call-with-values

					(lambda () (new-domain->old-domain (+ lower-0 1) lower-1))

				      old-indexer)

				    base)

				 0))

		(increment-1 (if (< (+ lower-1 1) upper-1)

				 (- (call-with-values

					(lambda () (new-domain->old-domain lower-0 (+ lower-1 1)))

				      old-indexer)

				    base)

				 0)))

	   (build-indexer-2 base lower-0 lower-1 increment-0 increment-1)))

    ((3) (let* ((lower-0 (##Interval-lower-bound new-domain 0))

		(lower-1 (##Interval-lower-bound new-domain 1))

		(lower-2 (##Interval-lower-bound new-domain 2))

		(upper-0 (##Interval-upper-bound new-domain 0))

		(upper-1 (##Interval-upper-bound new-domain 1))

		(upper-2 (##Interval-upper-bound new-domain 2))

		(base (call-with-values

			  (lambda () (new-domain->old-domain lower-0 lower-1 lower-2))

			old-indexer))

		(increment-0 (if (< (+ lower-0 1) upper-0)

				 (- (call-with-values

					(lambda () (new-domain->old-domain (+ lower-0 1) lower-1 lower-2))

				      old-indexer)

				    base)

				 0))

		(increment-1 (if (< (+ lower-1 1) upper-1)

				 (- (call-with-values

					(lambda () (new-domain->old-domain lower-0 (+ lower-1 1) lower-2))

				      old-indexer)

				    base)

				 0))

		(increment-2 (if (< (+ lower-2 1) upper-2)

				 (- (call-with-values

					(lambda () (new-domain->old-domain lower-0 lower-1 (+ lower-2 1)))

				      old-indexer)

				    base)

				 0)))

	   (build-indexer-3 base lower-0 lower-1 lower-2 increment-0 increment-1 increment-2)))

    ((4) (let* ((lower-0 (##Interval-lower-bound new-domain 0))

		(lower-1 (##Interval-lower-bound new-domain 1))

		(lower-2 (##Interval-lower-bound new-domain 2))

		(lower-3 (##Interval-lower-bound new-domain 3))

		(upper-0 (##Interval-upper-bound new-domain 0))

		(upper-1 (##Interval-upper-bound new-domain 1))

		(upper-2 (##Interval-upper-bound new-domain 2))

		(upper-3 (##Interval-upper-bound new-domain 3))

		(base (call-with-values

			  (lambda () (new-domain->old-domain lower-0 lower-1 lower-2 lower-3))

			old-indexer))

		(increment-0 (if (< (+ lower-0 1) upper-0)

				 (- (call-with-values

					(lambda () (new-domain->old-domain (+ lower-0 1) lower-1 lower-2 lower-3))

				      old-indexer)

				    base)

				 0))

		(increment-1 (if (< (+ lower-1 1) upper-1)

				 (- (call-with-values

					(lambda () (new-domain->old-domain lower-0 (+ lower-1 1) lower-2 lower-3))

				      old-indexer)

				    base)

				 0))

		(increment-2 (if (< (+ lower-2 1) upper-2)

				 (- (call-with-values

					(lambda () (new-domain->old-domain lower-0 lower-1 (+ lower-2 1) lower-3))

				      old-indexer)

				    base)

				 0))

		(increment-3 (if (< (+ lower-3 1) upper-3)

				 (- (call-with-values

					(lambda () (new-domain->old-domain lower-0 lower-1 lower-2 (+ lower-3 1)))

				      old-indexer)

				    base)

				 0)))

	   (build-indexer-4 base lower-0 lower-1 lower-2 lower-3 increment-0 increment-1 increment-2 increment-3)))

    (else

     (let* ((lower-bounds (##Interval-lower-bounds->list new-domain))

	    (upper-bounds (##Interval-upper-bounds->list new-domain))

	    (base (call-with-values

		      (lambda () (apply new-domain->old-domain lower-bounds))

		    old-indexer))

	    (increments (let ((increments   (map (lambda (x) 0) lower-bounds))

			      (lower-bounds (map (lambda (x) x) lower-bounds)))

			  (let loop ((l lower-bounds)

				     (u upper-bounds)

				     (i increments)

				     (base base))

			    (if (null? l)

				increments

				(let ((new-base

				       (if (< (+ (car l) 1)

					      (car u))

					   (begin

					     (set-car! l (+ (car l) 1))

					     (let ((new-base (call-with-values

								 (lambda () (apply new-domain->old-domain lower-bounds))

							       old-indexer)))

					       (set-car! i (- new-base base))

					       new-base))

					   base)))

				  (loop (cdr l)

					(cdr u)

					(cdr i)

					new-base)))))))

       (build-indexer-generic base lower-bounds increments)))))

#|

In the next function old-indexer is a 1-1 mapping from old-domain to [0,N) for some N.