compound.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
   3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   4 ;;; Since 1991, ANSI was finally finished.  Edited for ANSI Common Lisp.
   5
   6 ;;; compound -- Compound data and element-wise mapping functions
   7 ;;;
   8 ;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   9 ;;; unrestricted use.
  10 ;;;
  11
  12 ;;;
  13 ;;; Package Setup
  14 ;;;
  15
  16 (in-package :cl-user)
  17
  18 (defpackage :lisp-stat-compound-data
  19   (:use :common-lisp
  20         :lisp-stat-object-system
  21         :lisp-stat-types)
  22   (:import-from :lisp-stat-fastmap fastmap)
  23   (:shadowing-import-from :lisp-stat-object-system
  24                           slot-value
  25                           call-next-method call-method)
  26   (:export compound-data-p *compound-data-proto*
  27            compound-object-p
  28            compound-data-seq compound-data-length
  29
  30            element-list element-seq
  31
  32            sort-data order rank
  33
  34            recursive-map-elements map-elements
  35
  36            repeat
  37            ;; export matrix-related functionality (not sure??)
  38
  39            check-sequence
  40            get-next-element make-next-element set-next-element
  41            sequencep iseq
  42
  43            ;; maybe?
  44            ordered-nneg-seq
  45            select
  46
  47            which
  48            ;; vector differences
  49            difference rseq
  50            ))
  51
  52 (in-package :lisp-stat-compound-data)
  53
  54 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  55 ;;;
  56 ;;;                    Internal Support Functions
  57 ;;;
  58 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  59
  60 (defun cmpndp (x)
  61   "Predicate to determine if argument is compound. Most common
  62 non-compound types are checked first."
  63   (declare (inline numberp symbolp stringp consp arrayp array-total-size))
  64   (cond ((or (numberp x) (symbolp x) (stringp x)) nil)
  65          ((or (consp x) (and (arrayp x) (< 0 (array-total-size x)))) t)
  66          (t (compound-object-p x))))
  67
  68 (defun find-compound-data (list)
  69   "Returns first compound data item in LIST or NIL if there is none."
  70   (dolist (x list) (if (cmpndp x) (return x))))
  71
  72 (defun any-compound-elements (seq)
  73   "Checks for a compound element."
  74   (cond ((consp seq) (dolist (x seq) (if (cmpndp x) (return x))))
  75         ((vectorp seq)
  76          (let ((n (length seq)))
  77            (declare (fixnum n))
  78            (dotimes (i n)
  79              (declare (fixnum i))
  80              (let ((x (aref seq i)))
  81                (if (cmpndp x) (return x))))))
  82         (t (error "argument must be a list or vector"))))
  83
  84 (defun compound-data-sequence (x)
  85   "Returns sequence of data values for X."
  86   (declare (inline consp vectorp arrayp make-array array-total-size))
  87   (cond
  88    ((or (consp x) (vectorp x)) x)
  89    ((arrayp x) (make-array (array-total-size x) :displaced-to x))
  90    (t (send x :data-seq))))
  91
  92 (defmacro sequence-type (x) `(if (consp ,x) 'list 'vector))
  93
  94 (defun make-compound-data (shape sequence)
  95 "Construct a compound data item to match the shape of the first
  96 argument."
  97   (let ((n (length (compound-data-sequence shape))))
  98     (if (/= n (length sequence)) (error "compound data not the same shape"))
  99     (cond
 100      ((consp shape) (if (consp sequence) sequence (coerce sequence 'list)))
 101      ((vectorp shape)
 102       (if (vectorp sequence) sequence (coerce sequence 'vector)))
 103      ((arrayp shape)
 104       (make-array (array-dimensions shape)
 105                   :displaced-to (coerce sequence 'vector)))
 106      (t (send shape :make-data sequence)))))
 107
 108 (defun make-circle (x)
 109   "Make a circular list of one element."
 110   (declare (inline cons rplacd))
 111   (let ((x (cons x nil)))
 112     (rplacd x x)
 113     x))
 114
 115 (defun check-compound (x)
 116   "Signals an error if X is not compound."
 117   (if (not (cmpndp x)) (error "not a compound data item - ~a" x)))
 118
 119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 120 ;;;
 121 ;;; MAP-ELEMENTS function
 122 ;;; Applies a function to arguments. If all arguments are simple (i. e.
 123 ;;; not compound) then MAP-ELEMENTS acts like funcall. Otherwise all
 124 ;;; compound arguments must be of the same shape and simple arguments
 125 ;;; are treated as if they were compound arguments of the appropriate
 126 ;;; shape. This is implemented by replacin all simple arguments by
 127 ;;; circular lists of one element.
 128 ;;;
 129 ;;; This implementation uses FASTMAP, a version of MAP that is assumed
 130 ;;; to
 131 ;;;
 132 ;;;     a) work reasonable fast on any combination of lists and vectors
 133 ;;;        as its arguments
 134 ;;;
 135 ;;;     b) not hang if at least one of its arguments is not a circular
 136 ;;;        list.
 137 ;;;
 138 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 139
 140 (defun fixup-map-elements-arglist (args)
 141   (do* ((args args (rest args))
 142         (x (car args) (car args)))
 143        ((null args))
 144     (declare (inline car))
 145     (setf (car args)
 146           (if (cmpndp x) (compound-data-sequence x) (make-circle x)))))
 147
 148 (defun map-elements (fcn &rest args)
 149 "Args: (fcn &rest args)
 150 Applies FCN elementwise. If no arguments are compound MAP-ELEMENTS
 151 acts like FUNCALL. Compound arguments must all be the same shape. Non
 152 compound arguments, in the presence of compound ones, are treated as
 153 if they were of the same shape as the compound items with constant data
 154 values."
 155   (let ((first-compound (find-compound-data args)))
 156     (cond ((null first-compound) (apply fcn args))
 157           (t (fixup-map-elements-arglist args)
 158              (let* ((seq (compound-data-sequence first-compound))
 159                     (type (sequence-type seq)))
 160                (make-compound-data first-compound
 161                                    (apply #'fastmap type fcn args)))))))
 162
 163 (defun recursive-map-elements (base-fcn fcn &rest args)
 164 "Args: (base-fcn fcn &rest args)
 165 The same idea as MAP-ELEMENTS, except arguments are in a list and the
 166 base and recursive cases can use different functions. Modified to check
 167 for second level of compounding and use base-fcn if there is none."
 168   (let ((first-compound (find-compound-data args)))
 169     (cond ((null first-compound) (apply base-fcn args))
 170           (t (fixup-map-elements-arglist args)
 171              (let* ((seq (compound-data-sequence first-compound))
 172                     (type (sequence-type seq))
 173                     (f (if (any-compound-elements seq) fcn base-fcn)))
 174                (make-compound-data first-compound
 175                                    (apply #'fastmap type f args)))))))
 176
 177
 178 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 179 ;;;;
 180 ;;;;             Public Predicate and Accessor Functions
 181 ;;;;
 182 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 183
 184 (defun compound-data-p (x)
 185 "Args: (x)
 186 Returns T if X is a compound data item, NIL otherwise."
 187   (cmpndp x))
 188
 189 (defun compound-data-seq (x)
 190 "Args (x)
 191 Returns data sequence in X."
 192   (check-compound x)
 193   (compound-data-sequence x))
 194
 195 (defun compound-data-length (x)
 196 "Args (x)
 197 Returns length of data sequence in X."
 198   (check-compound x)
 199   (length (compound-data-sequence x)))
 200
 201 (defun compound-data-shape (x)
 202   "Needed but undefined??"
 203   x)
 204
 205
 206 (defun element-list (x)
 207   (cond
 208    ((compound-data-p x)
 209     (let ((x (concatenate 'list (compound-data-seq x)))) ; copies sequence
 210       (cond
 211        ((any-compound-elements x)
 212         (do ((next x (rest next)))
 213             ((not (consp next)))
 214           (setf (first next) (element-list (first next))))
 215         (do ((result (first x))
 216              (last (last (first x)))
 217              (next (rest x) (rest next)))
 218             ((not (consp next)) result)
 219           (setf (rest last) (first next))
 220           (setf last (last (first next)))))
 221        (t x))))
 222    (t (list x))))
 223
 224 (defun element-seq (x)
 225 "Args: (x)
 226 Returns sequence of the elements of compound item X."
 227   (check-compound x)
 228   (let ((seq (compound-data-seq x)))
 229     (if (any-compound-elements seq) (element-list seq) seq)))
 230
 231 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 232 ;;;;
 233 ;;;;                     Compound Data Objects
 234 ;;;;
 235 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 236
 237 (defproto *compound-data-proto*)
 238
 239 (defmeth *compound-data-proto* :data-length (&rest args) nil)
 240 (defmeth *compound-data-proto* :data-seq (&rest args) nil)
 241 (defmeth *compound-data-proto* :make-data (&rest args) nil)
 242 (defmeth *compound-data-proto* :select-data (&rest args) nil)
 243
 244 (defun compound-object-p (x) (kind-of-p x *compound-data-proto*))
 245
 246
 247
 248 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 249 ;;;;
 250 ;;;;                         Sorting Functions
 251 ;;;;
 252 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 253
 254 (defun sort-data (x)
 255 "Args: (sequence)
 256 Returns a sequence with the numbers or strings in the sequence X in order."
 257   (flet ((less (x y) (if (numberp x) (< x y) (string-lessp x y))))
 258     (stable-sort (copy-seq (compound-data-seq x)) #'less)))
 259
 260 (defun order (x)
 261 "Args (x)
 262 Returns a sequence of the indices of elements in the sequence of numbers
 263 or strings X in order."
 264   (let* ((seq (compound-data-seq x))
 265          (type (if (consp seq) 'list 'vector))
 266          (i -1))
 267     (flet ((entry (x) (setf i (+ i 1)) (list x i))
 268            (less (a b)
 269                  (let ((x (first a))
 270                        (y (first b)))
 271                    (if (numberp x) (< x y) (string-lessp x y)))))
 272       (let ((sorted-seq (stable-sort (map type #'entry seq) #'less)))
 273         (map type #'second sorted-seq)))))
 274
 275 ;; this isn't destructive -- do we document destructive only, or any
 276 ;; variant?
 277 (defun rank (x)
 278 "Args (x)
 279 Returns a sequence with the elements of the list or array of numbers or
 280 strings X replaced by their ranks."
 281   (let ((ranked-seq (order (order x))))
 282     (make-compound-data
 283      ;; compound-data-shape is undefined?
 284      (compound-data-shape x) ranked-seq)))
 285
 286
 287
 288 ;;;
 289 ;;; REPEAT function
 290 ;;;
 291
 292 (defun repeat (a b)
 293 "Args: (vals times)
 294 Repeats VALS. If TIMES is a number and VALS is a non-null, non-array atom,
 295 a list of length TIMES with all elements eq to VALS is returned. If VALS
 296 is a list and TIMES is a number then VALS is appended TIMES times. If
 297 TIMES is a list of numbers then VALS must be a list of equal length and
 298 the simpler version of repeat is mapped down the two lists.
 299 Examples: (repeat 2 5)                 returns (2 2 2 2 2)
 300           (repeat '(1 2) 3)            returns (1 2 1 2 1 2)
 301           (repeat '(4 5 6) '(1 2 3))   returns (4 5 5 6 6 6)
 302           (repeat '((4) (5 6)) '(2 3)) returns (4 4 5 6 5 6 5 6)"
 303   (cond ((compound-data-p b)
 304          (let* ((reps (coerce (compound-data-seq (map-elements #'repeat a b))
 305                               'list))
 306                 (result (first reps))
 307                 (tail (last (first reps))))
 308            (dolist (next (rest reps) result)
 309                    (when next
 310                          (setf (rest tail) next)
 311                          (setf tail (last next))))))
 312         (t (let* ((a (if (compound-data-p a)
 313                          (coerce (compound-data-seq a) 'list)
 314                          (list a)))
 315                   (result nil))
 316              (dotimes (i b result)
 317                       (let ((next (copy-list a)))
 318                         (if result (setf (rest (last next)) result))
 319                         (setf result next)))))))
 320 ;;;
 321 ;;; WHICH function
 322 ;;;
 323
 324 (defun which (x)
 325 "Args: (x)
 326 Returns a list of the indices where elements of sequence X are not NIL."
 327   (let ((x (list (compound-data-seq x)))
 328         (result nil)
 329         (tail nil))
 330     (flet ((add-result (x)
 331              (if result (setf (rest tail) (list x)) (setf result (list x)))
 332              (setf tail (if tail (rest tail) result)))
 333            (get-next-element (seq-list i)
 334              (cond ((consp (first seq-list))
 335                     (let ((elem (first (first seq-list))))
 336                       (setf (first seq-list) (rest (first seq-list)))
 337                       elem))
 338                    (t (aref (first seq-list) i)))))
 339           (let ((n (length (first x))))
 340             (dotimes (i n result)
 341                      (if (get-next-element x i) (add-result i)))))))
 342
 343
 344 ;;; Sequences are part of ANSI CL, being a supertype of vector and
 345 ;;; list (ordered set of things).
 346 ;;;
 347 ;;; Need to use the interenal structure when possible -- silly to be
 348 ;;; redundant!  However, this means we need to understand what
 349 ;;; sequences were intending to do, which I'm not clear on yet.
 350
 351 ;;; The original ordering, object-wise, was to have compound
 352 ;;; functionality passed into sequences, into other data sources.
 353 ;;; However, at this point, we will see about inverting this and
 354 ;;; having basic data types pushed through compound, to simplify
 355 ;;; packaging.
 356
 357 ;;;                      Type Checking Functions
 358
 359 (defun check-sequence (a)
 360   ;; FIXME:AJR: does this handle consp as well?  (Luke had an "or"
 361   ;; with consp).
 362   (if (not (typep a 'sequence))
 363       (error "not a sequence - ~s" a)))
 364
 365 ;;;                       Sequence Element Access
 366
 367
 368 ;;; (elt x i) -- NOT.  This is more like "pop".
 369 (defun get-next-element (x i)
 370   "Get element i from seq x.  FIXME: not really??"
 371   (let ((myseq (first x)))
 372     (if (consp myseq)
 373         (let ((elem (first myseq)))
 374           (setf (first x) (rest myseq))
 375           elem)
 376       (aref myseq i))))
 377
 378 ;;; (setf (elt x i) v)
 379 (defun set-next-element (x i v)
 380   (let ((seq (first x)))
 381     (cond ((consp seq)
 382            (setf (first seq) v)
 383            (setf (first x) (rest seq)))
 384           (t (setf (aref seq i) v)))))
 385
 386 (defun make-next-element (x) (list x))
 387
 388
 389 ;;;                         Sequence Functions
 390
 391
 392 ;; to prevent breakage.
 393 (defmacro sequencep (x)
 394   (typep x 'sequence))
 395
 396 (defun iseq (a &optional b)
 397 "Args: (n &optional m)
 398 Generate a sequence of consecutive integers from a to b.
 399 With one argumant returns a list of consecutive integers from 0 to N - 1.
 400 With two returns a list of consecutive integers from N to M.
 401 Examples: (iseq 4) returns (0 1 2 3)
 402           (iseq 3 7)  returns (3 4 5 6 7)
 403           (iseq 3 -3) returns (3 2 1 0 -1 -2 -3)"
 404   (if b
 405       (let ((n (+ 1 (abs (- b a))))
 406             (x nil))
 407         (dotimes (i n x)
 408                  (setq x (cons (if (< a b) (- b i) (+ b i)) x))))
 409       (cond
 410        ((= 0 a) nil)
 411        ((< a 0) (iseq (+ a 1) 0))
 412        ((< 0 a) (iseq 0 (- a 1))))))
 413
 414 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 415 ;;;;
 416 ;;;;               Subset Selection and Mutation Functions
 417 ;;;;
 418 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 419
 420
 421 (defun old-rowmajor-index (index indices dim olddim)
 422   "translate row major index in resulting subarray to row major index
 423 in the original array."
 424   (declare (fixnum index))
 425   (let ((rank (length dim))
 426         (face 1)
 427         (oldface 1)
 428         (oldindex 0))
 429     (declare (fixnum rank face oldface))
 430
 431     (dotimes (i rank)
 432       (declare (fixnum i))
 433       (setf face (* face (aref dim i)))
 434       (setf oldface (* oldface (aref olddim i))))
 435
 436     (dotimes (i rank)
 437       (declare (fixnum i))
 438       (setf face (/ face (aref dim i)))
 439       (setf oldface (/ oldface (aref olddim i)))
 440       (incf oldindex
 441             (* oldface (aref (aref indices i) (floor (/ index face))))) ;;*** is this floor really needed???
 442       (setf index (rem index face)))
 443     oldindex))
 444
 445 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 446 ;;;;
 447 ;;;;               Subset Selection and Mutation Functions
 448 ;;;;
 449 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 450
 451 (defun subarray-select (a indexlist &optional (values nil set_values))
 452   "extract or set subarray for the indices from a displaced array."
 453   (let ((indices nil)
 454         (index)
 455         (dim)
 456         (vdim)
 457         (data)
 458         (result_data)
 459         (olddim)
 460         (result)
 461         (rank 0)
 462         (n 0)
 463         (k 0))
 464     (declare (fixnum rank n))
 465
 466     (if (or (sequencep a) (not (arrayp a))) (error "not an array - ~a" a))
 467     (if (not (listp indexlist))  (error "bad index list - ~a" indices))
 468     (if (/= (length indexlist)  (array-rank a))
 469         (error "wrong number of indices"))
 470
 471     (setf indices (coerce indexlist 'vector))
 472
 473     (setf olddim (coerce (array-dimensions a) 'vector))
 474
 475     ;; compute the result dimension vector and fix up the indices
 476     (setf rank (array-rank a))
 477     (setf dim (make-array rank))
 478     (dotimes (i rank)
 479       (declare (fixnum i))
 480       (setf index (aref indices i))
 481       (setf n (aref olddim i))
 482       (setf index (if (fixnump index) (vector index) (coerce index 'vector)))
 483       (setf k (length index))
 484       (dotimes (j k)
 485         (declare (fixnum j))
 486         (if (<= n (check-nonneg-fixnum (aref index j)))
 487           (error "index out of bounds - ~a" (aref index j)))
 488         (setf (aref indices i) index))
 489       (setf (aref dim i) (length index)))
 490
 491     ;; set up the result or check the values
 492     (let ((dim-list (coerce dim 'list)))
 493       (cond
 494        (set_values
 495         (cond
 496          ((compound-data-p values)
 497           (if (or (not (arrayp values)) (/= rank (array-rank values)))
 498             (error "bad values array - ~a" values))
 499           (setf vdim (coerce (array-dimensions values) 'vector))
 500           (dotimes (i rank)
 501             (declare (fixnum i))
 502             (if (/= (aref vdim i) (aref dim i))
 503               (error "bad value array dimensions - ~a" values)))
 504           (setf result values))
 505          (t (setf result (make-array dim-list :initial-element values)))))
 506        (t (setf result (make-array dim-list)))))
 507
 508     ;; compute the result or set the values
 509     (setf data (compound-data-seq a))
 510     (setf result_data (compound-data-seq result))
 511     (setf n (length result_data))
 512     (dotimes (i n)
 513       (declare (fixnum i))
 514       (setf k (old-rowmajor-index i indices dim olddim))
 515       (if (or (> 0 k) (>= k (length data))) (error "index out of range"))
 516       (if set_values
 517         (setf (aref data k) (aref result_data i))
 518         (setf (aref result_data i) (aref data k))))
 519
 520     result))
 521
 522
 523 ;;;; is x an ordered sequence of nonnegative positive integers?
 524 (defun ordered-nneg-seq(x)
 525   ;; FIXME -- sbcl warning about unreachable code, might be a logic error here.
 526   (if (sequencep x)
 527       (let ((n (length x))
 528             (cx (make-next-element x))
 529             (m 0))
 530         (dotimes (i n t)
 531           (let ((elem (check-nonneg-fixnum (get-next-element cx i))))
 532             (if (> m elem) (return nil) (setf m elem)))))))
 533
 534 ;;;; select or set the subsequence corresponding to the specified indices
 535 (defun sequence-select(x indices &optional (values nil set-values))
 536   ;; FIXME -- sbcl warning about unreachable code, might be a logic error here.
 537   (let ((rlen 0)
 538         (dlen 0)
 539         (vlen 0)
 540         (data nil)
 541         (result nil))
 542     (declare (fixnum rlen dlen vlen))
 543
 544     ;; Check the input data
 545     (check-sequence x)
 546     (check-sequence indices)
 547     (if set-values (check-sequence values))
 548
 549     ;; Find the data sizes
 550     (setf data (if (ordered-nneg-seq indices) x (coerce x 'vector)))
 551     (setf dlen (length data))
 552     (setf rlen (length indices))
 553     (when set-values
 554       (setf vlen (length values))
 555       (if (/= vlen rlen) (error "value and index sequences do not match")))
 556
 557     ;; set up the result/value sequence
 558     (setf result
 559           (if set-values
 560               values
 561               (make-sequence (if (listp x) 'list 'vector) rlen)))
 562
 563     ;; get or set the sequence elements
 564     (if set-values
 565       (do ((nextx x)
 566            (cr (make-next-element result))
 567            (ci (make-next-element indices))
 568            (i 0 (+ i 1))
 569            (j 0)
 570            (index 0))
 571           ((>= i rlen))
 572         (declare (fixnum i j index))
 573         (setf index (get-next-element ci i))
 574         (if (<= dlen index) (error "index out of range - ~a" index))
 575         (let ((elem (get-next-element cr i)))
 576           (cond
 577            ((listp x)
 578             (when (> j index)
 579               (setf j 0)
 580               (setf nextx x))
 581             (do ()
 582                 ((not (and (< j index) (consp nextx))))
 583               (incf j 1)
 584               (setf nextx (rest nextx)))
 585             (setf (first nextx) elem))
 586            (t (setf (aref x index) elem)))))
 587       (do ((nextx data)
 588            (cr (make-next-element result))
 589            (ci (make-next-element indices))
 590            (i 0 (+ i 1))
 591            (j 0)
 592            (index 0)
 593            (elem nil))
 594           ((>= i rlen))
 595         (declare (fixnum i j index))
 596         (setf index (get-next-element ci i))
 597         (if (<= dlen index) (error "index out of range - ~a" index))
 598         (cond
 599          ((listp data) ;; indices must be ordered
 600           (do ()
 601               ((not (and (< j index) (consp nextx))))
 602             (incf j 1)
 603             (setf nextx (rest nextx)))
 604           (setf elem (first nextx)))
 605          (t (setf elem (aref data index))))
 606         (set-next-element cr i elem)))
 607
 608     result))
 609
 610 ;;;
 611 ;;; SELECT function
 612 ;;;
 613
 614 (defun select (x &rest args)
 615 "Args: (a &rest indices)
 616 A can be a list or an array. If A is a list and INDICES is a single number
 617 then the appropriate element of A is returned. If  is a list and INDICES is
 618 a list of numbers then the sublist of the corresponding elements is returned.
 619 If A in an array then the number of INDICES must match the ARRAY-RANK of A.
 620 If each index is a number then the appropriate array element is returned.
 621 Otherwise the INDICES must all be lists of numbers and the corresponding
 622 submatrix of A is returned. SELECT can be used in setf."
 623   (cond
 624    ((every #'fixnump args)
 625     (if (listp x) (nth (first args) x) (apply #'aref x args)))
 626    ((sequencep x) (sequence-select x (first args)))
 627    (t (subarray-select x args))))
 628
 629
 630 ;; Built in SET-SELECT (SETF method for SELECT)
 631 (defun set-select (x &rest args)
 632   (let ((indices (butlast args))
 633         (values (first (last args))))
 634     (cond
 635      ((sequencep x)
 636       (if (not (consp indices)) (error "bad indices - ~a" indices))
 637       (let* ((indices (first indices))
 638              (i-list (if (fixnump indices) (list indices) indices))
 639              (v-list (if (fixnump indices) (list values) values)))
 640         (sequence-select x i-list v-list)))
 641      ((arrayp x)
 642       (subarray-select x indices values))
 643      (t (error "bad argument type - ~a" x)))
 644     values))
 645
 646 (defsetf select set-select)
 647
 648 ;;;;
 649 ;;;; Basic Sequence Operations
 650 ;;;;
 651
 652 (defun difference (x)
 653 "Args: (x)
 654 Returns differences for a sequence X."
 655   (let ((n (length x)))
 656     (- (select x (iseq 1 (1- n))) (select x (iseq 0 (- n 2))))))
 657
 658 (defun rseq (a b num)
 659 "Args: (a b num)
 660 Returns a list of NUM equally spaced points starting at A and ending at B."
 661   (+ a (* (values-list (iseq 0 (1- num))) (/ (float (- b a)) (1- num)))))
 662
 663