Language/OS - Multiplatform Resource Library/LANGUAGE OS.iso/scheme/pcscheme/geneva/pcscheme.exe/SAMPLES/ARRAY.S

Internet Message Format | 1993-05-13 | 11.4 KB

From news.unige.ch!scsing.switch.ch!ira.uka.de!yale.edu!qt.cs.utexas.edu!cs.utexas.edu!sun-barr!olivea!mintaka.lcs.mit.edu!mintaka.lcs.mit.edu!alan Mon Feb 22 10:27:37 1993
Path: news.unige.ch!scsing.switch.ch!ira.uka.de!yale.edu!qt.cs.utexas.edu!cs.utexas.edu!sun-barr!olivea!mintaka.lcs.mit.edu!mintaka.lcs.mit.edu!alan
From: [e-mail] (Alan Bawden)
Newsgroups: comp.lang.scheme
Subject: Re: multi-dimensional arrays in R5RS?
Message-ID: <[e-mail]>
Date: 21 Feb 93 05:56:39 GMT
References: <[e-mail]>
    <[e-mail]>
    <[e-mail]>
    <[e-mail]>
    <[e-mail]>
Sender: [e-mail]
Organization: ITS Preservation Society
Lines: 326
In-Reply-To: [e-mail]'s message of 20 Feb 93 06:09:14 GMT

In article <[e-mail]>
[e-mail] (david carlton) writes:

   In article <[e-mail]>,
   "John Lacey" <[e-mail]> 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                               [e-mail]
Work: MIT Room NE43-510, 545 Tech. Sq., Cambridge, MA 02139     [phone]
Home: 29 Reed St., Cambridge, MA 02140                 [phone]