diff --git a/lib/srfi/179.scm b/lib/srfi/179.scm index 0b60192d..b9495ab0 100644 --- a/lib/srfi/179.scm +++ b/lib/srfi/179.scm @@ -1,194 +1,7 @@ -;; Miscellaneous Functions -(define (translation? x) - (and (vector? x) (not (vector-empty? x)) (vector-every exact-integer? x))) - -(define (permutation? x) - (and (translation? x) - (let* ((len (vector-length x)) - (seen (make-u1vector len))) - (let lp ((i 0)) - (or (>= i len) - (and (< -1 (vector-ref x i) len) - (zero? (u1vector-ref seen (vector-ref x i))) - (begin - (u1vector-set! seen (vector-ref x i) 1) - (lp (+ i 1))))))))) - -(define (same-dimensions? ls) - (or (null? ls) - (null? (cdr ls)) - (and (equal? (array-dimension (car ls)) (array-dimension (cadr ls))) - (same-dimensions? (cdr ls))))) - -;; Intervals - -(define-record-type Interval - (%%make-interval lb ub) - interval? - (lb interval-lb) - (ub interval-ub)) - -(define (%make-interval lo hi) - (assert (and (translation? lo) - (translation? hi) - (= (vector-length lo) (vector-length hi)) - (vector-every < lo hi))) - (%%make-interval lo hi)) - -(define (make-interval x . o) - (if (pair? o) - (%make-interval x (car o)) - (%make-interval (make-vector (vector-length x) 0) x))) - -(define (interval-dimension iv) - (vector-length (interval-lb iv))) - -(define (interval-lower-bound iv i) (vector-ref (interval-lb iv) i)) -(define (interval-upper-bound iv i) (vector-ref (interval-ub iv) i)) -(define (interval-lower-bounds->list iv) (vector->list (interval-lb iv))) -(define (interval-upper-bounds->list iv) (vector->list (interval-ub iv))) -(define (interval-lower-bounds->vector iv) (vector-copy (interval-lb iv))) -(define (interval-upper-bounds->vector iv) (vector-copy (interval-ub iv))) - -(define (interval= iv1 iv2) - (assert (and (interval? iv1) (interval? iv2))) - (equal? iv1 iv2)) - -(define (interval-volume iv) - (vector-fold (lambda (acc lower upper) (* acc (- upper lower))) - 1 - (interval-lb iv) (interval-ub iv))) - -(define (interval-subset? iv1 iv2) - (assert (and (interval? iv1) (interval? iv2) - (= (interval-dimension iv1) (interval-dimension iv2)))) - (and (vector-every >= (interval-lb iv1) (interval-lb iv2)) - (vector-every <= (interval-ub iv1) (interval-ub iv2)))) - -(define (interval-contains-multi-index? iv i0 . o) - (assert (interval? iv)) - (let ((i (list->vector (cons i0 o)))) - (assert (and (= (interval-dimension iv) (vector-length i)) - (vector-every integer? i))) - (and (vector-every >= i (interval-lb iv)) - (vector-every < i (interval-ub iv))))) - -(define (interval-projections iv rd) - (values (make-interval (vector-copy (interval-lb iv) 0 rd) - (vector-copy (interval-ub iv) 0 rd)) - (make-interval (vector-copy (interval-lb iv) rd) - (vector-copy (interval-ub iv) rd)))) - -(define (rev-index-next! rev-index rev-lowers rev-uppers) - (cond - ((null? rev-index) #f) - ((< (caar rev-index) (- (car rev-uppers) 1)) - (set-car! (car rev-index) (+ 1 (caar rev-index))) - #t) - (else - (set-car! (car rev-index) (car rev-lowers)) - (rev-index-next! (cdr rev-index) (cdr rev-lowers) (cdr rev-uppers))))) - -(define (interval-fold kons knil iv) - (case (interval-dimension iv) - ((1) - (let ((end (interval-upper-bound iv 0))) - (do ((i (interval-lower-bound iv 0) (+ i 1)) - (acc knil (kons acc i))) - ((>= i end) acc)))) - ((2) - (let ((end0 (interval-upper-bound iv 0)) - (start1 (interval-lower-bound iv 1)) - (end1 (interval-upper-bound iv 1))) - (do ((i (interval-lower-bound iv 0) (+ i 1)) - (acc knil - (do ((j start1 (+ j 1)) - (acc acc (kons acc i j))) - ((>= j end1) acc)))) - ((>= i end0) acc)))) - (else - (let* ((rev-lowers (reverse (interval-lower-bounds->list iv))) - (rev-uppers (reverse (interval-upper-bounds->list iv))) - (multi-index (interval-lower-bounds->list iv)) - (rev-index (pair-fold cons '() multi-index))) - (let lp ((acc knil)) - (let ((acc (apply kons acc multi-index))) - (if (rev-index-next! rev-index rev-lowers rev-uppers) - (lp acc) - acc))))))) - -(define (interval-for-each f iv) - (interval-fold (lambda (acc . multi-index) (apply f multi-index)) #f iv) - (if #f #f)) - -(define (interval-dilate iv lower-diffs upper-diffs) - (assert (= (interval-dimension iv) - (vector-length lower-diffs) - (vector-length upper-diffs))) - (make-interval (vector-map + (interval-lb iv) lower-diffs) - (vector-map + (interval-ub iv) upper-diffs))) - -(define (interval-intersect iv0 . o) - (let ((ls (cons iv0 o))) - (assert (and (every interval? ls) - (or (null? o) (apply = (map interval-dimension ls))))) - (let ((lower (apply vector-map max (map interval-lb ls))) - (upper (apply vector-map min (map interval-ub ls)))) - (and (vector-every < lower upper) - (make-interval lower upper))))) - -(define (interval-translate iv translation) - (assert (translation? translation)) - (interval-dilate iv translation translation)) - -(define (interval-permute iv perm) - (assert (and (interval? iv) (permutation? perm))) - (let* ((len (interval-dimension iv)) - (lower (make-vector len)) - (upper (make-vector len))) - (assert (= len (vector-length perm))) - (do ((i 0 (+ i 1))) - ((>= i len) (make-interval lower upper)) - (vector-set! lower i (interval-lower-bound iv (vector-ref perm i))) - (vector-set! upper i (interval-upper-bound iv (vector-ref perm i)))))) - -(define (interval-rotate iv dim) - (let ((lower (interval-lb iv)) - (upper (interval-ub iv))) - (make-interval (vector-append (vector-copy lower dim) - (vector-copy lower 0 dim)) - (vector-append (vector-copy upper dim) - (vector-copy upper 0 dim))))) - -(define (interval-scale iv scales) - (assert (and (interval? iv) - (vector? scales) - (= (interval-dimension iv) (vector-length scales)) - (vector-every exact-integer? scales) - (vector-every positive? scales))) - (make-interval - (vector-map (lambda (u s) (exact (ceiling (/ u s)))) - (interval-ub iv) - scales))) - -(define (interval-cartesian-product iv0 . o) - (make-interval (apply vector-append (map interval-lb (cons iv0 o))) - (apply vector-append (map interval-ub (cons iv0 o))))) - -;; Storage Classes - -(define-record-type Storage-Class - (make-storage-class getter setter checker maker copier length default) - storage-class? - (getter storage-class-getter) - (setter storage-class-setter) - (checker storage-class-checker) - (maker storage-class-maker) - (copier storage-class-copier) - (length storage-class-length) - (default storage-class-default)) +;; Homogeneous storage classes +;; Define a storage class with an optimized -copy! (define-syntax define-storage-class (syntax-rules () ((define-storage-class name ref set elt? make len default) @@ -206,11 +19,6 @@ (set to i (ref from j)))))) len default))))) -(define generic-storage-class - (make-storage-class - vector-ref vector-set! (lambda (x) #t) make-vector - vector-copy! vector-length #f)) - (define-storage-class s8-storage-class s8vector-ref s8vector-set! s8? make-s8vector s8vector-length 0) @@ -250,260 +58,10 @@ (define-storage-class c128-storage-class c128vector-ref c128vector-set! c128? make-c128vector c128vector-length 0) +;; TODO: implement (define f8-storage-class #f) (define f16-storage-class #f) -;; Arrays - -(define-record-type Array - (%%make-array domain getter setter storage body coeffs indexer safe? adjacent?) - array? - (domain array-domain) - (getter array-getter) - (setter array-setter %array-setter-set!) - (storage array-storage-class) - (body array-body) - (coeffs array-coeffs) - (indexer array-indexer) - (safe? array-safe?) - (adjacent? array-adjacent? array-adjacent?-set!)) - -(define (%make-array domain getter setter storage body coeffs - indexer safe? adjacent?) - (assert (and (interval? domain) - (procedure? getter) - (or (not setter) (procedure? setter)) - (or (not storage) (storage-class? storage)))) - (%%make-array - domain getter setter storage body coeffs indexer safe? adjacent?)) - -(define (make-array domain getter . o) - (assert (and (interval? domain) (procedure? getter))) - (%make-array domain getter (and (pair? o) (car o)) #f #f #f #f #f #f)) - -(define (array-dimension a) - (interval-dimension (array-domain a))) - -(define (mutable-array? x) - (and (array? x) (array-setter x) #t)) - -(define specialized-array-default-safe? - (make-parameter #f (lambda (x) (assert (boolean? x)) x))) - -(define specialized-array-default-mutable? - (make-parameter #t (lambda (x) (assert (boolean? x)) x))) - -(define (array-ref array . multi-index) - (apply (array-getter array) multi-index)) - -(define (array-set! array val . multi-index) - (apply (array-setter array) val multi-index)) - -(define (specialized-getter body indexer getter) - (lambda multi-index - (getter body (apply indexer multi-index)))) - -(define (specialized-setter body indexer setter) - (lambda (val . multi-index) - (setter body (apply indexer multi-index) val))) - -;; Indexing - -(define (indexer->coeffs indexer domain . o) - (let* ((verify? (and (pair? o) (car o))) - (res (make-vector (+ 1 (interval-dimension domain)) 0)) - (multi-index (interval-lower-bounds->list domain)) - (base (apply indexer multi-index))) - (vector-set! res 0 base) - (let lp ((i 1) - (ls multi-index) - (offset base) - (count 0)) - (cond - ((null? ls) - (if (and verify? (zero? count)) - (lp 1 multi-index offset (+ count 1)) - res)) - ((= (+ 1 (interval-lower-bound domain (- i 1))) - (interval-upper-bound domain (- i 1))) - (lp (+ i 1) (cdr ls) offset count)) - (else - (let ((dir (if (and (> count 0) - (= (+ (car ls) 1) - (interval-upper-bound domain (- i 1)))) - -1 - 1))) - (set-car! ls (+ (car ls) dir)) - (let* ((offset2 (apply indexer multi-index)) - (coeff (* dir (- offset2 offset)))) - (cond - ((> count 0) - (and (= coeff (vector-ref res i)) - (lp (+ i 1) (cdr ls) offset2 count))) - (else - (vector-set! res i coeff) - (vector-set! res 0 (- (vector-ref res 0) - (* coeff - (interval-lower-bound domain (- i 1))))) - (lp (+ i 1) (cdr ls) offset2 count)))))))))) - -(define (coeffs->indexer coeffs domain) - (case (vector-length coeffs) - ((2) - (let ((a (vector-ref coeffs 0)) - (b (vector-ref coeffs 1))) - (lambda (x) (+ a (* b x))))) - ((3) - (let ((a (vector-ref coeffs 0)) - (b (vector-ref coeffs 1)) - (c (vector-ref coeffs 2))) - (lambda (x y) (+ a (* b x) (* c y))))) - ((4) - (let ((a (vector-ref coeffs 0)) - (b (vector-ref coeffs 1)) - (c (vector-ref coeffs 2)) - (d (vector-ref coeffs 3))) - (lambda (x y z) (+ a (* b x) (* c y) (* d z))))) - (else - (lambda multi-index - (let ((lim (vector-length coeffs))) - (let lp ((ls multi-index) - (i 1) - (res (vector-ref coeffs 0))) - (cond - ((null? ls) - (if (< i lim) - (error "multi-index too short for domain" multi-index domain) - res)) - ((>= i lim) - (error "multi-index too long for domain" multi-index domain)) - (else - (lp (cdr ls) - (+ i 1) - (+ res (* (car ls) (vector-ref coeffs i)))))))))))) - -(define (default-coeffs domain) - (let* ((dim (interval-dimension domain)) - (res (make-vector (+ 1 dim)))) - (vector-set! res 0 0) - (vector-set! res dim 1) - (let lp ((i (- dim 1)) - (scale 1)) - (cond - ((< i 0) - res) - ((= (+ 1 (interval-lower-bound domain i)) - (interval-upper-bound domain i)) - (vector-set! res (+ i 1) 0) - (lp (- i 1) scale)) - (else - (let ((coeff (* scale (- (interval-upper-bound domain i) - (interval-lower-bound domain i))))) - (vector-set! res (+ i 1) scale) - (vector-set! res 0 (- (vector-ref res 0) - (* scale (interval-lower-bound domain i)))) - (lp (- i 1) coeff))))))) - -(define (default-indexer domain) - (coeffs->indexer (default-coeffs domain) domain)) - -;; Converts the raw integer index to the multi-index in domain that -;; would map to it using the default indexer (i.e. iterating over the -;; possible multi-indices in domain in lexicographic order would -;; produce 0 through volume-1). -(define (invert-default-index domain raw-index) - (let lp ((index raw-index) - (i 0) - (scale (/ (interval-volume domain) - (max 1 - (- (interval-upper-bound domain 0) - (interval-lower-bound domain 0))))) - (res '())) - (cond - ((>= (+ i 1) (interval-dimension domain)) - (reverse (cons (+ index (interval-lower-bound domain i)) res))) - (else - (let ((digit (quotient index scale))) - (lp (- index (* digit scale)) - (+ i 1) - (/ scale - (max 1 - (- (interval-upper-bound domain (+ i 1)) - (interval-lower-bound domain (+ i 1))))) - (cons (+ digit - (interval-lower-bound domain i)) - res))))))) - -;; Specialized arrays - -(define (%make-specialized domain storage body coeffs indexer - safe? mutable? adjacent?) - (%make-array - domain - (specialized-getter body indexer (storage-class-getter storage)) - (and mutable? - (specialized-setter body indexer (storage-class-setter storage))) - storage - body - coeffs - indexer - safe? - adjacent?)) - -(define (make-specialized-array domain . o) - (let* ((storage (if (pair? o) (car o) generic-storage-class)) - (safe? (if (and (pair? o) (pair? (cdr o))) - (cadr o) - (specialized-array-default-safe?))) - (body ((storage-class-maker storage) - (interval-volume domain) - (storage-class-default storage))) - (coeffs (default-coeffs domain)) - (indexer (coeffs->indexer coeffs domain))) - (assert (boolean? safe?)) - (%make-specialized domain storage body coeffs indexer safe? #t #t))) - -(define (specialized-array? x) - (and (array? x) (array-storage-class x) #t)) - -(define (compute-array-elements-in-order? array) - (let ((indexer (array-indexer array))) - (call-with-current-continuation - (lambda (return) - (interval-fold - (lambda (prev . multi-index) - (let ((i (apply indexer multi-index))) - (if (and prev (not (= i (+ prev 1)))) - (return #f) - i))) - #f - (array-domain array)) - #t)))) - -(define (array-elements-in-order? array) - (assert (specialized-array? array)) - (let ((res (array-adjacent? array))) - (when (eq? res 'unknown) - (set! res (compute-array-elements-in-order? array)) - (array-adjacent?-set! array res)) - res)) - -(define (specialized-array-share array new-domain project) - (assert (and (specialized-array? array) (interval? new-domain))) - (let* ((body (array-body array)) - (coeffs - (indexer->coeffs - (lambda multi-index - (call-with-values - (lambda () (apply project multi-index)) - (array-indexer array))) - new-domain)) - (indexer - (coeffs->indexer coeffs new-domain)) - (storage (array-storage-class array))) - (%make-specialized new-domain storage body coeffs indexer - (array-safe? array) (array-setter array) 'unknown))) - ;; Array transformations (define (array-copy array . o) @@ -727,6 +285,12 @@ (op (apply getter1 (take multi-index dim1)) (apply getter2 (drop multi-index dim1))))))) +(define (same-dimensions? ls) + (or (null? ls) + (null? (cdr ls)) + (and (equal? (array-dimension (car ls)) (array-dimension (cadr ls))) + (same-dimensions? (cdr ls))))) + (define (array-map f array . arrays) (make-array (array-domain array) (let* ((ls (cons array arrays)) @@ -815,7 +379,7 @@ (let ((getter (array-getter source)) (setter (array-setter destination))) (cond - ((equal? (array-domain destination) (array-domain source)) + ((interval= (array-domain destination) (array-domain source)) (interval-for-each (case (array-dimension destination) ((1) (lambda (i) (setter (getter i) i))) @@ -829,15 +393,12 @@ (assert (and (array-elements-in-order? destination) (equal? (interval-volume (array-domain destination)) (interval-volume (array-domain source))))) - (let* ((dst-domain (array-domain destination)) - (rev-lowers (reverse (interval-lower-bounds->list dst-domain))) - (rev-uppers (reverse (interval-upper-bounds->list dst-domain))) - (dst-index (list-copy (interval-lower-bounds->list dst-domain))) - (rev-index (pair-fold cons '() dst-index))) + (let* ((ivc (interval-cursor (array-domain destination))) + (dst-index (interval-cursor-get ivc))) (interval-for-each (lambda multi-index (apply setter (apply getter multi-index) dst-index) - (rev-index-next! rev-index rev-lowers rev-uppers)) + (interval-cursor-next! ivc)) (array-domain source))))) destination)) diff --git a/lib/srfi/179.sld b/lib/srfi/179.sld index cd95060d..ccc07d94 100644 --- a/lib/srfi/179.sld +++ b/lib/srfi/179.sld @@ -4,6 +4,7 @@ (scheme vector) (scheme sort) (srfi 160 base) + (srfi 179 base) (chibi assert)) (export ;; Miscellaneous Functions diff --git a/lib/srfi/179/base.scm b/lib/srfi/179/base.scm new file mode 100644 index 00000000..4442530a --- /dev/null +++ b/lib/srfi/179/base.scm @@ -0,0 +1,469 @@ +;; Miscellaneous Functions + +(define (translation? x) + (and (vector? x) (not (vector-empty? x)) (vector-every exact-integer? x))) + +(define (permutation? x) + (and (translation? x) + (let* ((len (vector-length x)) + (seen (make-vector len 0))) + (let lp ((i 0)) + (or (>= i len) + (and (< -1 (vector-ref x i) len) + (zero? (vector-ref seen (vector-ref x i))) + (begin + (vector-set! seen (vector-ref x i) 1) + (lp (+ i 1))))))))) + +;; Intervals + +(define-record-type Interval + (%%make-interval lb ub) + interval? + (lb interval-lb) + (ub interval-ub)) + +(define (%make-interval lo hi) + (assert (and (translation? lo) + (translation? hi) + (= (vector-length lo) (vector-length hi)) + (vector-every < lo hi))) + (%%make-interval lo hi)) + +(define (make-interval x . o) + (if (pair? o) + (%make-interval x (car o)) + (%make-interval (make-vector (vector-length x) 0) x))) + +(define (interval-dimension iv) + (vector-length (interval-lb iv))) + +(define (interval-lower-bound iv i) (vector-ref (interval-lb iv) i)) +(define (interval-upper-bound iv i) (vector-ref (interval-ub iv) i)) +(define (interval-lower-bounds->list iv) (vector->list (interval-lb iv))) +(define (interval-upper-bounds->list iv) (vector->list (interval-ub iv))) +(define (interval-lower-bounds->vector iv) (vector-copy (interval-lb iv))) +(define (interval-upper-bounds->vector iv) (vector-copy (interval-ub iv))) + +(define (interval= iv1 iv2) + (assert (and (interval? iv1) (interval? iv2))) + (and (equal? (interval-lb iv1) (interval-lb iv2)) + (equal? (interval-ub iv1) (interval-ub iv2)))) + +(define (interval-volume iv) + (vector-fold (lambda (acc lower upper) (* acc (- upper lower))) + 1 + (interval-lb iv) (interval-ub iv))) + +(define (interval-subset? iv1 iv2) + (assert (and (interval? iv1) (interval? iv2) + (= (interval-dimension iv1) (interval-dimension iv2)))) + (and (vector-every >= (interval-lb iv1) (interval-lb iv2)) + (vector-every <= (interval-ub iv1) (interval-ub iv2)))) + +(define (interval-contains-multi-index? iv i0 . o) + (assert (interval? iv)) + (let ((i (list->vector (cons i0 o)))) + (assert (and (= (interval-dimension iv) (vector-length i)) + (vector-every integer? i))) + (and (vector-every >= i (interval-lb iv)) + (vector-every < i (interval-ub iv))))) + +(define (interval-projections iv rd) + (values (make-interval (vector-copy (interval-lb iv) 0 rd) + (vector-copy (interval-ub iv) 0 rd)) + (make-interval (vector-copy (interval-lb iv) rd) + (vector-copy (interval-ub iv) rd)))) + +(define (rev-index-next! rev-index rev-lowers rev-uppers) + (cond + ((null? rev-index) #f) + ((< (caar rev-index) (- (car rev-uppers) 1)) + (set-car! (car rev-index) (+ 1 (caar rev-index))) + #t) + (else + (set-car! (car rev-index) (car rev-lowers)) + (rev-index-next! (cdr rev-index) (cdr rev-lowers) (cdr rev-uppers))))) + +(define (interval-cursor iv) + (let* ((rev-lowers (reverse (interval-lower-bounds->list iv))) + (rev-uppers (reverse (interval-upper-bounds->list iv))) + (multi-index (interval-lower-bounds->list iv)) + (rev-index (pair-fold cons '() multi-index))) + (vector multi-index rev-index rev-lowers rev-uppers))) + +(define (interval-cursor-get ivc) + (vector-ref ivc 0)) + +(define (interval-cursor-next! ivc) + (and (rev-index-next! (vector-ref ivc 1) + (vector-ref ivc 2) + (vector-ref ivc 3)) + (vector-ref ivc 0))) + +(define (interval-cursor-next ivc) + (let* ((multi-index (list-copy (vector-ref ivc 0))) + (ivc (vector multi-index + (pair-fold cons '() multi-index) + (vector-ref ivc 2) + (vector-ref ivc 3)))) + (and (rev-index-next! (vector-ref ivc 1) + (vector-ref ivc 2) + (vector-ref ivc 3)) + (values ivc (vector-ref ivc 0))))) + +(define (interval-fold kons knil iv) + (case (interval-dimension iv) + ((1) + (let ((end (interval-upper-bound iv 0))) + (do ((i (interval-lower-bound iv 0) (+ i 1)) + (acc knil (kons acc i))) + ((>= i end) acc)))) + ((2) + (let ((end0 (interval-upper-bound iv 0)) + (start1 (interval-lower-bound iv 1)) + (end1 (interval-upper-bound iv 1))) + (do ((i (interval-lower-bound iv 0) (+ i 1)) + (acc knil + (do ((j start1 (+ j 1)) + (acc acc (kons acc i j))) + ((>= j end1) acc)))) + ((>= i end0) acc)))) + (else + (let ((ivc (interval-cursor iv))) + (let lp ((acc knil)) + (let ((acc (apply kons acc (interval-cursor-get ivc)))) + (if (interval-cursor-next! ivc) + (lp acc) + acc))))))) + +(define (interval-for-each f iv) + (interval-fold (lambda (acc . multi-index) (apply f multi-index)) #f iv) + (if #f #f)) + +(define (interval-dilate iv lower-diffs upper-diffs) + (assert (= (interval-dimension iv) + (vector-length lower-diffs) + (vector-length upper-diffs))) + (make-interval (vector-map + (interval-lb iv) lower-diffs) + (vector-map + (interval-ub iv) upper-diffs))) + +(define (interval-intersect iv0 . o) + (let ((ls (cons iv0 o))) + (assert (and (every interval? ls) + (or (null? o) (apply = (map interval-dimension ls))))) + (let ((lower (apply vector-map max (map interval-lb ls))) + (upper (apply vector-map min (map interval-ub ls)))) + (and (vector-every < lower upper) + (make-interval lower upper))))) + +(define (interval-translate iv translation) + (assert (translation? translation)) + (interval-dilate iv translation translation)) + +(define (interval-permute iv perm) + (assert (and (interval? iv) (permutation? perm))) + (let* ((len (interval-dimension iv)) + (lower (make-vector len)) + (upper (make-vector len))) + (assert (= len (vector-length perm))) + (do ((i 0 (+ i 1))) + ((>= i len) (make-interval lower upper)) + (vector-set! lower i (interval-lower-bound iv (vector-ref perm i))) + (vector-set! upper i (interval-upper-bound iv (vector-ref perm i)))))) + +(define (interval-rotate iv dim) + (let ((lower (interval-lb iv)) + (upper (interval-ub iv))) + (make-interval (vector-append (vector-copy lower dim) + (vector-copy lower 0 dim)) + (vector-append (vector-copy upper dim) + (vector-copy upper 0 dim))))) + +(define (interval-scale iv scales) + (assert (and (interval? iv) + (vector? scales) + (= (interval-dimension iv) (vector-length scales)) + (vector-every exact-integer? scales) + (vector-every positive? scales))) + (make-interval + (vector-map (lambda (u s) (exact (ceiling (/ u s)))) + (interval-ub iv) + scales))) + +(define (interval-cartesian-product iv0 . o) + (make-interval (apply vector-append (map interval-lb (cons iv0 o))) + (apply vector-append (map interval-ub (cons iv0 o))))) + +;; Storage Classes + +(define-record-type Storage-Class + (make-storage-class getter setter checker maker copier length default) + storage-class? + (getter storage-class-getter) + (setter storage-class-setter) + (checker storage-class-checker) + (maker storage-class-maker) + (copier storage-class-copier) + (length storage-class-length) + (default storage-class-default)) + +(define generic-storage-class + (make-storage-class + vector-ref vector-set! (lambda (x) #t) make-vector + vector-copy! vector-length #f)) + +;; Parameters + +;; Note safety is ignored in this implementation. +(define specialized-array-default-safe? + (make-parameter #f (lambda (x) (assert (boolean? x)) x))) + +(define specialized-array-default-mutable? + (make-parameter #t (lambda (x) (assert (boolean? x)) x))) + +;; Arrays + +(define-record-type Array + (%%make-array domain getter setter storage body coeffs indexer safe? adjacent?) + array? + (domain array-domain) + (getter array-getter) + (setter array-setter %array-setter-set!) + (storage array-storage-class) + (body array-body) + (coeffs array-coeffs) + (indexer array-indexer) + (safe? array-safe?) + (adjacent? array-adjacent? array-adjacent?-set!)) + +(define (%make-array domain getter setter storage body coeffs + indexer safe? adjacent?) + (assert (and (interval? domain) + (procedure? getter) + (or (not setter) (procedure? setter)) + (or (not storage) (storage-class? storage)))) + (%%make-array + domain getter setter storage body coeffs indexer safe? adjacent?)) + +(define (make-array domain getter . o) + (assert (and (interval? domain) (procedure? getter))) + (%make-array domain getter (and (pair? o) (car o)) #f #f #f #f #f #f)) + +(define (array-dimension a) + (interval-dimension (array-domain a))) + +(define (mutable-array? x) + (and (array? x) (array-setter x) #t)) + +(define (array-ref array . multi-index) + (apply (array-getter array) multi-index)) + +(define (array-set! array val . multi-index) + (apply (array-setter array) val multi-index)) + +(define (specialized-getter body indexer getter) + (lambda multi-index + (getter body (apply indexer multi-index)))) + +(define (specialized-setter body indexer setter) + (lambda (val . multi-index) + (setter body (apply indexer multi-index) val))) + + +;; Indexing + +(define (indexer->coeffs indexer domain . o) + (let* ((verify? (and (pair? o) (car o))) + (res (make-vector (+ 1 (interval-dimension domain)) 0)) + (multi-index (interval-lower-bounds->list domain)) + (base (apply indexer multi-index))) + (vector-set! res 0 base) + (let lp ((i 1) + (ls multi-index) + (offset base) + (count 0)) + (cond + ((null? ls) + (if (and verify? (zero? count)) + (lp 1 multi-index offset (+ count 1)) + res)) + ((= (+ 1 (interval-lower-bound domain (- i 1))) + (interval-upper-bound domain (- i 1))) + (lp (+ i 1) (cdr ls) offset count)) + (else + (let ((dir (if (and (> count 0) + (= (+ (car ls) 1) + (interval-upper-bound domain (- i 1)))) + -1 + 1))) + (set-car! ls (+ (car ls) dir)) + (let* ((offset2 (apply indexer multi-index)) + (coeff (* dir (- offset2 offset)))) + (cond + ((> count 0) + (and (= coeff (vector-ref res i)) + (lp (+ i 1) (cdr ls) offset2 count))) + (else + (vector-set! res i coeff) + (vector-set! res 0 (- (vector-ref res 0) + (* coeff + (interval-lower-bound domain (- i 1))))) + (lp (+ i 1) (cdr ls) offset2 count)))))))))) + +(define (coeffs->indexer coeffs domain) + (case (vector-length coeffs) + ((2) + (let ((a (vector-ref coeffs 0)) + (b (vector-ref coeffs 1))) + (lambda (x) (+ a (* b x))))) + ((3) + (let ((a (vector-ref coeffs 0)) + (b (vector-ref coeffs 1)) + (c (vector-ref coeffs 2))) + (lambda (x y) (+ a (* b x) (* c y))))) + ((4) + (let ((a (vector-ref coeffs 0)) + (b (vector-ref coeffs 1)) + (c (vector-ref coeffs 2)) + (d (vector-ref coeffs 3))) + (lambda (x y z) (+ a (* b x) (* c y) (* d z))))) + (else + (lambda multi-index + (let ((lim (vector-length coeffs))) + (let lp ((ls multi-index) + (i 1) + (res (vector-ref coeffs 0))) + (cond + ((null? ls) + (if (< i lim) + (error "multi-index too short for domain" multi-index domain) + res)) + ((>= i lim) + (error "multi-index too long for domain" multi-index domain)) + (else + (lp (cdr ls) + (+ i 1) + (+ res (* (car ls) (vector-ref coeffs i)))))))))))) + +(define (default-coeffs domain) + (let* ((dim (interval-dimension domain)) + (res (make-vector (+ 1 dim)))) + (vector-set! res 0 0) + (vector-set! res dim 1) + (let lp ((i (- dim 1)) + (scale 1)) + (cond + ((< i 0) + res) + ((= (+ 1 (interval-lower-bound domain i)) + (interval-upper-bound domain i)) + (vector-set! res (+ i 1) 0) + (lp (- i 1) scale)) + (else + (let ((coeff (* scale (- (interval-upper-bound domain i) + (interval-lower-bound domain i))))) + (vector-set! res (+ i 1) scale) + (vector-set! res 0 (- (vector-ref res 0) + (* scale (interval-lower-bound domain i)))) + (lp (- i 1) coeff))))))) + +(define (default-indexer domain) + (coeffs->indexer (default-coeffs domain) domain)) + +;; Converts the raw integer index to the multi-index in domain that +;; would map to it using the default indexer (i.e. iterating over the +;; possible multi-indices in domain in lexicographic order would +;; produce 0 through volume-1). +(define (invert-default-index domain raw-index) + (let lp ((index raw-index) + (i 0) + (scale (/ (interval-volume domain) + (max 1 + (- (interval-upper-bound domain 0) + (interval-lower-bound domain 0))))) + (res '())) + (cond + ((>= (+ i 1) (interval-dimension domain)) + (reverse (cons (+ index (interval-lower-bound domain i)) res))) + (else + (let ((digit (quotient index scale))) + (lp (- index (* digit scale)) + (+ i 1) + (/ scale + (max 1 + (- (interval-upper-bound domain (+ i 1)) + (interval-lower-bound domain (+ i 1))))) + (cons (+ digit + (interval-lower-bound domain i)) + res))))))) + +;; Specialized arrays + +(define (%make-specialized domain storage body coeffs indexer + safe? mutable? adjacent?) + (%make-array + domain + (specialized-getter body indexer (storage-class-getter storage)) + (and mutable? + (specialized-setter body indexer (storage-class-setter storage))) + storage + body + coeffs + indexer + safe? + adjacent?)) + +(define (make-specialized-array domain . o) + (let* ((storage (if (pair? o) (car o) generic-storage-class)) + (safe? (if (and (pair? o) (pair? (cdr o))) + (cadr o) + (specialized-array-default-safe?))) + (body ((storage-class-maker storage) + (interval-volume domain) + (storage-class-default storage))) + (coeffs (default-coeffs domain)) + (indexer (coeffs->indexer coeffs domain))) + (assert (boolean? safe?)) + (%make-specialized domain storage body coeffs indexer safe? #t #t))) + +(define (specialized-array? x) + (and (array? x) (array-storage-class x) #t)) + +(define (compute-array-elements-in-order? array) + (let ((indexer (array-indexer array))) + (call-with-current-continuation + (lambda (return) + (interval-fold + (lambda (prev . multi-index) + (let ((i (apply indexer multi-index))) + (if (and prev (not (= i (+ prev 1)))) + (return #f) + i))) + #f + (array-domain array)) + #t)))) + +(define (array-elements-in-order? array) + (assert (specialized-array? array)) + (let ((res (array-adjacent? array))) + (when (eq? res 'unknown) + (set! res (compute-array-elements-in-order? array)) + (array-adjacent?-set! array res)) + res)) + +(define (specialized-array-share array new-domain project) + (assert (and (specialized-array? array) (interval? new-domain))) + (let* ((body (array-body array)) + (coeffs + (indexer->coeffs + (lambda multi-index + (call-with-values + (lambda () (apply project multi-index)) + (array-indexer array))) + new-domain)) + (indexer + (coeffs->indexer coeffs new-domain)) + (storage (array-storage-class array))) + (%make-specialized new-domain storage body coeffs indexer + (array-safe? array) (array-setter array) 'unknown))) diff --git a/lib/srfi/179/base.sld b/lib/srfi/179/base.sld new file mode 100644 index 00000000..8f50178a --- /dev/null +++ b/lib/srfi/179/base.sld @@ -0,0 +1,41 @@ + +;;> The base array definitions of SRFI 179, plus some extra internal +;;> bindings. + +(define-library (srfi 179 base) + (import (scheme base) + (scheme list) + (scheme vector) + (chibi assert)) + (export + ;; Miscellaneous Functions + translation? permutation? + ;; Intervals + make-interval interval? interval-dimension interval-lb interval-ub + interval-lower-bound interval-upper-bound interval-lower-bounds->list + interval-upper-bounds->list interval-lower-bounds->vector + interval-upper-bounds->vector interval= interval-volume + interval-subset? interval-contains-multi-index? interval-projections + interval-for-each interval-dilate interval-intersect + interval-translate interval-permute interval-rotate + interval-scale interval-cartesian-product + ;; Indexing + indexer->coeffs coeffs->indexer default-indexer default-coeffs + invert-default-index interval-cursor interval-cursor-next! + interval-cursor-next interval-cursor-get interval-fold + ;; Storage Classes + make-storage-class storage-class? storage-class-getter + storage-class-setter storage-class-checker storage-class-maker + storage-class-copier storage-class-length storage-class-default + generic-storage-class + ;; Arrays + make-array array? array-domain array-getter array-dimension + mutable-array? array-setter specialized-array-default-safe? + specialized-array-default-mutable? make-specialized-array + specialized-array? array-storage-class array-indexer array-body + array-safe? array-adjacent? array-elements-in-order? + specialized-array-share array-ref array-set! + %make-specialized %array-setter-set! + specialized-getter specialized-setter + ) + (include "base.scm"))