diff options
Diffstat (limited to 'array.scm')
| -rw-r--r-- | array.scm | 279 |
1 files changed, 279 insertions, 0 deletions
diff --git a/array.scm b/array.scm new file mode 100644 index 0000000..3eecb7a --- /dev/null +++ b/array.scm @@ -0,0 +1,279 @@ +;;;;"array.scm" Arrays for Scheme +; Copyright (C) 1993 Alan Bawden +; +; Permission to copy this software, to redistribute it, 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. Users of this software agree to make their best efforts (a) to +; return to me any improvements or extensions that they make, so that +; these may be included in future releases; and (b) to inform me of +; noteworthy uses of this software. +; +; 3. I have made no warrantee 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. +; +; 4. 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. +; +; Alan Bawden +; MIT Room NE43-510 +; 545 Tech. Sq. +; Cambridge, MA 02139 +; Alan@LCS.MIT.EDU + +(require 'record) + +;(declare (usual-integrations)) + +(define array:rtd + (make-record-type "Array" + '(indexer ; Must be a -linear- function! + shape ; Inclusive bounds: ((lower upper) ...) + vector ; The actual contents + ))) + +(define array:indexer (record-accessor array:rtd 'indexer)) +(define array-shape (record-accessor array:rtd 'shape)) +(define array:vector (record-accessor array:rtd 'vector)) + +(define array? (record-predicate array:rtd)) + +(define (array-rank obj) + (if (array? obj) (length (array-shape obj)) 0)) + +(define (array-dimensions ra) + (map (lambda (ind) (if (zero? (car ind)) (cadr ind) ind)) + (array-shape ra))) + +(define array:construct + (record-constructor array:rtd '(shape vector indexer))) + +(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 (slib:error "array: Bad array dimension: " 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 shape + (make-vector size initial-value) + (array:optimize-linear-function indexer 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 new-shape + (array:vector array) + (array:optimize-linear-function + (lambda indices + (apply old-indexer (apply mapping indices))) + 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) + (slib:error "array: Bad indices for " 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)) + +(define (array-in-bounds? array . indices) + (array:in-bounds? array indices)) + +; Fast versions of ARRAY-REF and ARRAY-SET! that do no error checking, +; and don't cons intermediate lists of indices: + +(define (array-1d-ref a i0) + (vector-ref (array:vector a) ((array:indexer a) i0))) + +(define (array-2d-ref a i0 i1) + (vector-ref (array:vector a) ((array:indexer a) i0 i1))) + +(define (array-3d-ref a i0 i1 i2) + (vector-ref (array:vector a) ((array:indexer a) i0 i1 i2))) + +(define (array-1d-set! a v i0) + (vector-set! (array:vector a) ((array:indexer a) i0) v)) + +(define (array-2d-set! a v i0 i1) + (vector-set! (array:vector a) ((array:indexer a) i0 i1) v)) + +(define (array-3d-set! a v i0 i1 i2) + (vector-set! (array:vector a) ((array:indexer a) i0 i1 i2) v)) + +; 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 l) 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 l) + (let ((d (length l))) + (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 (map (lambda (x) 0) l)) + (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-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)))))) |