; From: alan@lcs.mit.edu (Alan Bawden)
; Newsgroups: comp.lang.scheme
; Subject: Re: multi-dimensional arrays in R5RS?
; Date: 21 Feb 93 05:56:39 GMT
; Organization: ITS Preservation Society
; In-Reply-To: carlton@husc11.harvard.edu's message of 20 Feb 93 06:09:14 GMT
; 
; In article <CARLTON.93Feb19220914@husc11.harvard.edu> 
; carlton@husc11.harvard.edu (david carlton) writes:
; 
;    In article <1993Feb19.142555.15696@news.cs.indiana.edu>,
;    "John Lacey" <johnl@cs.indiana.edu> writes:
;    > An implication of this approach to the issue at hand is that the
;    > printed representation would be a vector of vectors, but the
;    > compiler would see multidimensional array primitives, with all
;    > the optimization benefits of that.
;    This approach (which I think is a good one) says nothing about the
;    printed representation.  The library could provide functions to
;    manipulate multidimensional arrays without guaranteeing anything at
;    all about their printed representation or their representation.
; 
; Right.
; 
; The most peculiar thing to me about this discussion of arrays for Scheme is
; the seemingly widespread notion that arrays should be somehow identified
; with vectors that contain vectors.  (One common corollary of this belief
; seems to be that if you don't have native arrays, you -have- to use vectors
; of vectors -- which is inefficient.)
; 
; Well, let me point out that the difference between a vector and an array is
; really just some arithmetic performed on the indices.  If I'm writing a
; program that does a lot of manipulating of 4x4 matrices, I can do as well
; as any native array implementation by using 16 element vectors, and writing
; 
;   (vector-ref M (+ (* 4 i) j))
; 
; instead of
; 
;   (array-ref M i j)
; 
; If I write my loops so that, in fact, I process the elements of the matrix
; in the order that they are layed out in the vector, then a good compiler
; should be able to do just as well is it could if I used native arrays.
; 
; In order to drive this point home, here is an implementation of arrays for
; scheme that works by abstracting out the index arithmetic.  This makes it
; trivial to provide shared subarrays.  In fact, this approach gets you a lot
; more than you could ever get by simply using vectors of vectors.
; 
; (This code has -not- been thoroughly debugged -- near the end there are a
; bunch of little functions that I wrote by hand that have not been
; exhaustively tested and could easily contain typos.)
; 
------- Begin arrays.scm

; Arrays for Scheme
;
; Copyright (C) 1993 Alan Bawden
;
; You may use this code any way you like, as long as you don't charge money
; for it, remove this notice, or hold me liable for its results.

; The user interface consists of the following 5 functions
; 
;   (ARRAY? <object>)  =>  <boolean>
;   (MAKE-ARRAY <initial-value> <bound> <bound> ...)  =>  <array>
;   (ARRAY-REF <array> <index> <index> ...)  =>  <value>
;   (ARRAY-SET! <array> <new-value> <index> <index> ...)
;   (MAKE-SHARED-ARRAY <array> <mapper> <bound> <bound> ...)  =>  <array>
;
; When constructing an array, <bound> is either an inclusive range of
; indices expressed as a two element list, or an upper bound expressed as a
; single integer.  So
;
;   (make-array 'foo 3 3)
;
; and
;
;   (make-array 'foo '(0 2) '(0 2))
;
; are equivalent.
;
; MAKE-SHARED-ARRAY can be used to create shared subarrays of other arrays.
; The <mapper> is a function that translates coordinates in the new array
; into coordinates in the old array.  A <mapper> must be linear, and its
; range must stay within the bounds of the old array, but it can be
; otherwise arbitrary.  A simple example:
;
;   (define fred (make-array #F 8 8))
;   (define freds-diagonal
;     (make-shared-array fred (lambda (i) (list i i)) 8))
;   (array-set! freds-diagonal 'foo 3)
;   (array-ref fred 3 3)  =>  FOO
;   (define freds-center
;     (make-shared-array fred (lambda (i j) (list (+ 3 i) (+ 3 j))) 2 2))
;   (array-ref freds-center 0 0)  =>  FOO
;
; End of manual.

; If you comment out the bounds checking code, this is about as efficient
; as you could ask for without help from the compiler.
;
; An exercise left to the reader: implement the rest of APL.
;
; Functions not part of the user interface have names that start with
; "array/" (the poor man's package system).
;
; This assumes your Scheme has the usual record type package.

(define array/rtd
  (make-record-type "Array"
    '(vector
      indexer		; Must be a -linear- function!
      shape)))

(define array/vector (record-accessor array/rtd 'vector))
(define array/indexer (record-accessor array/rtd 'indexer))
(define array/shape (record-accessor array/rtd 'shape))

(define array? (record-predicate array/rtd))

(define array/construct
  (record-constructor array/rtd '(vector indexer shape)))

(define (array/compute-shape specs)
  (map (lambda (spec)
	 (cond ((and (integer? spec)
		     (< 0 spec))
		(list 0 (- spec 1)))
	       ((and (pair? spec)
		     (pair? (cdr spec))
		     (null? (cddr spec))
		     (integer? (car spec))
		     (integer? (cadr spec))
		     (<= (car spec) (cadr spec)))
		spec)
	       (else (error "Bad array dimension: ~S" spec))))
       specs))

(define (make-array initial-value . specs)
  (let ((shape (array/compute-shape specs)))
    (let loop ((size 1)
	       (indexer (lambda () 0))
	       (l (reverse shape)))
      (if (null? l)
	  (array/construct (make-vector size initial-value)
			   (array/optimize-linear-function indexer
							   (length shape))
			   shape)
	  (loop (* size (+ 1 (- (cadar l) (caar l))))
		(lambda (first-index . rest-of-indices)
		  (+ (* size (- first-index (caar l)))
		     (apply indexer rest-of-indices)))
		(cdr l))))))

(define (make-shared-array array mapping . specs)
  (let ((new-shape (array/compute-shape specs))
	(old-indexer (array/indexer array)))
    (let check ((indices '())
		(bounds (reverse new-shape)))
      (cond ((null? bounds)
	     (array/check-bounds array (apply mapping indices)))
	    (else
	     (check (cons (caar bounds) indices) (cdr bounds))
	     (check (cons (cadar bounds) indices) (cdr bounds)))))
    (array/construct (array/vector array)
		     (array/optimize-linear-function
		       (lambda indices
			 (apply old-indexer (apply mapping indices)))
		       (length new-shape))
		     new-shape)))

(define (array/in-bounds? array indices)
  (let loop ((indices indices)
	     (shape (array/shape array)))
    (if (null? indices)
	(null? shape)
	(let ((index (car indices)))
	  (and (not (null? shape))
	       (integer? index)
	       (<= (caar shape) index (cadar shape))
	       (loop (cdr indices) (cdr shape)))))))

(define (array/check-bounds array indices)
  (or (array/in-bounds? array indices)
      (error "Bad indices for ~S: ~S" array indices)))

(define (array-ref array . indices)
  (array/check-bounds array indices)
  (vector-ref (array/vector array)
	      (apply (array/indexer array) indices)))

(define (array-set! array new-value . indices)
  (array/check-bounds array indices)
  (vector-set! (array/vector array)
	       (apply (array/indexer array) indices)
	       new-value))

; STOP!  Do not read beyond this point on your first reading of
; this code -- you should simply assume that the rest of this file
; contains only the following single definition:
;
;   (define (array/optimize-linear-function f d) f)
;
; Of course everything would be pretty inefficient if this were really the
; case, but it isn't.  The following code takes advantage of the fact that
; you can learn everything there is to know from a linear function by
; simply probing around in its domain and observing its values -- then a
; more efficient equivalent can be constructed.

(define (array/optimize-linear-function f d)
  (cond ((= d 0)
	 (array/0d-c (f)))
	((= d 1)
	 (let ((c (f 0)))
	   (array/1d-c0 c (- (f 1) c))))
	((= d 2)
	 (let ((c (f 0 0)))
	   (array/2d-c01 c (- (f 1 0) c) (- (f 0 1) c))))
	((= d 3)
	 (let ((c (f 0 0 0)))
	   (array/3d-c012 c (- (f 1 0 0) c) (- (f 0 1 0) c) (- (f 0 0 1) c))))
	(else
	 (let* ((v (make-list d 0))
		(c (apply f v)))
	   (let loop ((p v)
		      (old-val c)
		      (coefs '()))
	     (cond ((null? p)
		    (array/Nd-c* c (reverse coefs)))
		   (else
		    (set-car! p 1)
		    (let ((new-val (apply f v)))
		      (loop (cdr p)
			    new-val
			    (cons (- new-val old-val) coefs))))))))))

; 0D cases:

(define (array/0d-c c)
  (lambda () c))

; 1D cases:

(define (array/1d-c c)
  (lambda (i0) (+ c i0)))

(define (array/1d-0 n0)
  (cond ((= 1 n0) +)
	(else (lambda (i0) (* n0 i0)))))

(define (array/1d-c0 c n0)
  (cond ((= 0 c) (array/1d-0 n0))
	((= 1 n0) (array/1d-c c))
	(else (lambda (i0) (+ c (* n0 i0))))))

; 2D cases:

(define (array/2d-0 n0)
  (lambda (i0 i1) (+ (* n0 i0) i1)))

(define (array/2d-1 n1)
  (lambda (i0 i1) (+ i0 (* n1 i1))))

(define (array/2d-c0 c n0)
  (lambda (i0 i1) (+ c (* n0 i0) i1)))

(define (array/2d-c1 c n1)
  (lambda (i0 i1) (+ c i0 (* n1 i1))))

(define (array/2d-01 n0 n1)
  (cond ((= 1 n0) (array/2d-1 n1))
	((= 1 n1) (array/2d-0 n0))
	(else (lambda (i0 i1) (+ (* n0 i0) (* n1 i1))))))

(define (array/2d-c01 c n0 n1)
  (cond ((= 0 c) (array/2d-01 n0 n1))
	((= 1 n0) (array/2d-c1 c n1))
	((= 1 n1) (array/2d-c0 c n0))
	(else (lambda (i0 i1) (+ c (* n0 i0) (* n1 i1))))))

; 3D cases:

(define (array/3d-01 n0 n1)
  (lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) i2)))

(define (array/3d-02 n0 n2)
  (lambda (i0 i1 i2) (+ (* n0 i0) i1 (* n2 i2))))

(define (array/3d-12 n1 n2)
  (lambda (i0 i1 i2) (+ i0 (* n1 i1) (* n2 i2))))

(define (array/3d-012 n0 n1 n2)
  (lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) (* n2 i2))))

(define (array/3d-c12 c n1 n2)
  (lambda (i0 i1 i2) (+ c i0 (* n1 i1) (* n2 i2))))

(define (array/3d-c02 c n0 n2)
  (lambda (i0 i1 i2) (+ c (* n0 i0) i1 (* n2 i2))))

(define (array/3d-c01 c n0 n1)
  (lambda (i0 i1 i2) (+ c (* n0 i0) (* n1 i1) i2)))

(define (array/3d-012 n0 n1 n2)
  (cond ((= 1 n0) (array/3d-12 n1 n2))
	((= 1 n1) (array/3d-02 n0 n2))
	((= 1 n2) (array/3d-01 n0 n1))
	(else (lambda (i0 i1 i2) (+ (* n0 i0) (* n1 i1) (* n2 i2))))))

(define (array/3d-c012 c n0 n1 n2)
  (cond ((= 0 c) (array/3d-012 n0 n1 n2))
	((= 1 n0) (array/3d-c12 c n1 n2))
	((= 1 n1) (array/3d-c02 c n0 n2))
	((= 1 n2) (array/3d-c01 c n0 n1))
	(else (lambda (i0 i1 i2) (+ c (* n0 i0) (* n1 i1) (* n2 i2))))))

; ND cases:

(define (array/Nd-* coefs)
  (lambda indices (apply + (map * coefs indices))))

(define (array/Nd-c* c coefs)
  (cond ((= 0 c) (array/Nd-* coefs))
	(else (lambda indices (apply + c (map * coefs indices))))))

------- End arrays.scm
--

--
Alan Bawden						       Alan@LCS.MIT.EDU
Work: MIT Room NE43-510, 545 Tech. Sq., Cambridge, MA 02139	 (617) 253-7328
Home: 29 Reed St., Cambridge, MA 02140				 (617) 492-7274