share/amatrix/amatrix.lisp

   1 ;; amatrix.lisp -- implement Maxima matrix via underlying Lisp array
   2 ;; copyright 2007 by Robert Dodier
   3 ;; I release this file under the terms of the GNU General Public License.
   4
   5 (setf (get '$amatrix '$present) t)
   6
   7 (declare-top (unspecial x))
   8
   9 (defun $amatrixmap (f &rest m)
  10   (let*
  11     ((nr ($nrows (car m)))
  12      (nc ($ncols (car m)))
  13      (m1 (mfuncall '$make_matrix nr nc))
  14      (a1 (get ($@-function m1 '$storage) 'storage_array)))
  15     (dotimes (i nr)
  16       (dotimes (j nc)
  17         (let*
  18           ((args
  19              (mapcar
  20                #'(lambda (m0)
  21                    (let
  22                      ((a0 (get ($@-function m0 '$storage) 'storage_array))
  23                       (ij (mfuncall '$compute_index0 m0 i j)))
  24                      (aref a0 ij)))
  25                m))
  26            (x (mapply1 f args f nil)))
  27           (setf (aref a1 (mfuncall '$compute_index0 m1 i j)) x))))
  28     m1))
  29
  30 (putprop '$amatrix 'amatrix-assign 'mset_extension_operator)
  31
  32 (defun amatrix-assign (e x)
  33   (let
  34     ((my-amatrix-name (caar e))
  35      (my-amatrix (symbol-value (caar e))))
  36     (cond
  37       ((= (length e) 2)
  38        (cond
  39          ((= ($nrows my-amatrix) 1)
  40           (amatrix-assign1 my-amatrix-name my-amatrix 1 (meval (cadr e)) x))
  41          ((= ($ncols my-amatrix) 1)
  42           (amatrix-assign1 my-amatrix-name my-amatrix (meval (cadr e)) 1 x))
  43          (t
  44            (merror "amatrix assignment: given one subscript, but expected two"))))
  45       ((= (length e) 3)
  46        (amatrix-assign1 (caar e) (symbol-value (caar e)) (meval (cadr e)) (meval (caddr e)) x))
  47       (t
  48         (merror "amatrix assignment: expected one or two subscripts")))))
  49
  50 (defun amatrix-assign1 (lhs aa i j x)
  51   (when (> (get ($@-function aa '$storage) 'refcount) 1)
  52     (let ((a (gensym)))
  53       (putprop a 1 'refcount)
  54       (putprop a (mfuncall '$copy_array (get ($@-function aa '$storage) 'storage_array)) 'storage_array)
  55       (decf (get ($@-function aa '$storage) 'refcount))
  56       (mrecord-assign `(($@) ,lhs $storage) a)))
  57   (cond
  58     ((and (integerp i) (integerp j))
  59      (amatrix-assign1-row-column aa i j x))
  60     ((and (eq i '$all) (integerp j))
  61      (amatrix-assign1-all-column aa j x))
  62     ((and (integerp i) (eq j '$all))
  63      (amatrix-assign1-row-all aa i x))
  64     ((and (eq i '$all) (eq j '$all))
  65      (amatrix-assign1-all-all aa x))
  66     ((and ($listp i) (integerp j))
  67      (amatrix-assign1-list-integer aa i j x))
  68     ((and (integerp i) ($listp j))
  69      (amatrix-assign1-integer-list aa i j x))
  70     ((and ($listp i) (eq j '$all))
  71      (amatrix-assign1-list-all aa i x))
  72     ((and (eq i '$all) ($listp j))
  73      (amatrix-assign1-all-list aa j x))
  74     ((and ($listp i) ($listp j))
  75      (amatrix-assign1-list-list aa i j x))
  76     (t
  77       `((mset) ((,aa array) ,i ,j) ,x))))
  78
  79 (defun amatrix-assign1-row-column (aa i j x)
  80   (let
  81     ((a (get ($@-function aa '$storage) 'storage_array))
  82      (ij (1- (mfuncall '$compute_index1 aa i j))))
  83     (setf (aref a ij) x)))
  84
  85 (defun amatrix-assign1-all-column (aa j x)
  86   (let ((m ($@-function aa '$nr)))
  87     (if ($amatrixp x)
  88       ;; MIGHT WANT TO ENSURE THAT X HAS SAME NUMBER OF ROWS AS AA AND EXACTLY ONE COLUMN
  89       (dotimes (i m)
  90         (amatrix-assign1-row-column aa (1+ i) j (mfuncall '$get_element x (1+ i) 1)))
  91       (dotimes (i m)
  92         (amatrix-assign1-row-column aa (1+ i) j x))))
  93   x)
  94
  95 (defun amatrix-assign1-row-all (aa i x)
  96   (let ((n ($@-function aa '$nc)))
  97     (if ($amatrixp x)
  98       ;; MIGHT WANT TO ENSURE THAT X HAS EXACTLY ONE ROW AND SAME NUMBER OF COLUMNS AS AA
  99       (dotimes (j n)
 100         (amatrix-assign1-row-column aa i (1+ j) (mfuncall '$get_element x 1 (1+ j))))
 101       (dotimes (j n)
 102         (amatrix-assign1-row-column aa i (1+ j) x))))
 103   x)
 104
 105 (defun amatrix-assign1-all-all (aa x)
 106   (let
 107     ((m ($@-function aa '$nr))
 108      (n ($@-function aa '$nc)))
 109     (if ($amatrixp x)
 110       ;; MIGHT WANT TO ENSURE THAT X HAS SAME NUMBER OF ROWS AND COLUMNS AS AA
 111       (dotimes (i m)
 112         (dotimes (j n)
 113           (amatrix-assign1-row-column aa (1+ i) (1+ j) (mfuncall '$get_element x (1+ i) (1+ j)))))
 114       (dotimes (i m)
 115         (dotimes (j n)
 116           (amatrix-assign1-row-column aa (1+ i) (1+ j) x)))))
 117   x)
 118
 119 (defun amatrix-assign1-list-integer (aa l j x)
 120   (amatrix-assign1-list-list aa l `((mlist) ,j) x))
 121
 122 (defun amatrix-assign1-integer-list (aa i l x)
 123   (amatrix-assign1-list-list aa `((mlist) ,i) l x))
 124
 125 (defun amatrix-assign1-list-all (aa l x)
 126   (let
 127     ((m ($length l))
 128      (n ($@-function aa '$nc)))
 129     (if ($amatrixp x)
 130       ;; MIGHT WANT TO ENSURE THAT X HAS SAME NUMBER OF ROWS AND COLUMNS AS AA
 131       (dotimes (i m)
 132         (dotimes (j n)
 133           (amatrix-assign1-row-column aa (nth (1+ i) l) (1+ j) (mfuncall '$get_element x (1+ i) (1+ j)))))
 134       (dotimes (i m)
 135         (dotimes (j n)
 136           (amatrix-assign1-row-column aa (nth (1+ i) l) (1+ j) x)))))
 137   x)
 138
 139 (defun amatrix-assign1-all-list (aa l x)
 140   (let
 141     ((m ($@-function aa '$nr))
 142      (n ($length l)))
 143     (if ($amatrixp x)
 144       ;; MIGHT WANT TO ENSURE THAT X HAS SAME NUMBER OF ROWS AND COLUMNS AS AA
 145       (dotimes (i m)
 146         (dotimes (j n)
 147           (amatrix-assign1-row-column aa (1+ i) (nth (1+ j) l) (mfuncall '$get_element x (1+ i) (1+ j)))))
 148       (dotimes (i m)
 149         (dotimes (j n)
 150           (amatrix-assign1-row-column aa (1+ i) (nth (1+ j) l) x)))))
 151   x)
 152
 153 (defun amatrix-assign1-list-list (aa l1 l2 x)
 154   (let
 155     ((m ($length l1))
 156      (n ($length l2)))
 157     (if ($amatrixp x)
 158       ;; MIGHT WANT TO ENSURE THAT X HAS SAME NUMBER OF ROWS AND COLUMNS AS AA
 159       (dotimes (i m)
 160         (dotimes (j n)
 161           (amatrix-assign1-row-column
 162             aa
 163             (nth (1+ i) l1)
 164             (nth (1+ j) l2)
 165             (mfuncall '$get_element x (1+ i) (1+ j)))))
 166       (dotimes (i m)
 167         (dotimes (j n)
 168           (amatrix-assign1-row-column
 169             aa
 170             (nth (1+ i) l1)
 171             (nth (1+ j) l2)
 172             x)))))
 173   x)
 174
 175 ;; UNFINISHED STUFF. SORRY, I'M WORKING ON IT.
 176
 177 (defun amatrix-assign1-amatrix-column (aa ai j x) (declare (ignore aa ai j x)))
 178
 179 (defun amatrix-assign1-row-amatrix (aa i aj x) (declare (ignore aa i aj x)))
 180
 181 (defun amatrix-assign1-amatrix-amatrix (aa ai aj x) (declare (ignore aa ai aj x)))
 182
 183 (defun amatrix-assign1-boolean-column (aa b j x) (declare (ignore aa b j x)))
 184
 185 (defun amatrix-assign1-row-boolean (aa i b x) (declare (ignore aa i b x)))
 186
 187 (defun amatrix-assign1-boolean-boolean (aa b1 b2 x) (declare (ignore aa b1 b2 x)))
 188
 189 (displa-def $amatrix dim-$amatrix)
 190
 191 ; CALL SIMPLIFYA IN DIM-$NEWMATRIX BECAUSE DIM-$MATRIX NEEDS TO SEE SIMP FLAG ... SIGH
 192
 193 (defun dim-$amatrix (form result)
 194   (let
 195     ((nr ($@-function form '$nr))
 196      (r0 ($@-function form '$r0))
 197      (rinc ($@-function form '$rinc))
 198      (nc ($@-function form '$nc))
 199      (c0 ($@-function form '$c0))
 200      (cinc ($@-function form '$cinc))
 201      (storage ($@-function form '$storage)))
 202     (if
 203       (and
 204         (integerp nr)
 205         (integerp r0)
 206         (integerp rinc)
 207         (integerp nc)
 208         (integerp c0)
 209         (integerp cinc)
 210         (symbolp storage)
 211         (arrayp (get storage 'storage_array)))
 212       (if (or (= nr 0) (= nc 0))
 213         (dimension-function '(($amatrix simp)) result)
 214         (dim-$matrix
 215           (simplifya
 216             ($genmatrix
 217               `((lambda) ((mlist) i j) (mfuncall '$get_element ,form i j))
 218               ($@-function form '$nr)
 219               ($@-function form '$nc))
 220             t)
 221           result))
 222       (dimension-function form result))))
 223
 224 (defun $amatrixp (a)
 225   (and (not ($atom a)) (eq ($op a) '$amatrix)))
 226
 227 (defun $nrows (a)
 228   (if ($amatrixp a)
 229     ($@-function a '$nr)
 230     `(($nrows) ,a)))
 231
 232 (defun $ncols (a)
 233   (if ($amatrixp a)
 234     ($@-function a '$nc)
 235     `(($ncols) ,a)))
 236
 237 (defun $amatrix_multiply (m1 m2)
 238   (let ((nc1 ($ncols m1)) (nr2 ($nrows m2)))
 239     (cond
 240       ((and (integerp nc1) (integerp nr2))
 241        (if (eql ($ncols m1) ($nrows m2))
 242          (let
 243            ((aa (gensym))
 244             (nn ($ncols m1))
 245             (nr ($nrows m1))
 246             (a1 (get ($@-function m1 '$storage) 'storage_array))
 247             (a1-inc ($@-function m1 '$cinc))
 248             (nc ($ncols m2))
 249             (a2 (get ($@-function m2 '$storage) 'storage_array))
 250             (a2-inc ($@-function m2 '$rinc)))
 251            (setf (get aa 'storage_array) (make-array (* nr nc) :initial-element 0))
 252            (dotimes (i nr)
 253              (dotimes (j nc)
 254                (let
 255                  ((a1-base (mfuncall '$compute_index0 m1 i 0))
 256                   (a2-base (mfuncall '$compute_index0 m2 0 j)))
 257                  (setf
 258                    (aref (get aa 'storage_array) (+ i (* j nr)))
 259                    (dot nn a1 a1-base a1-inc a2 a2-base a2-inc)))))
 260            `(($amatrix) ,nr 0 1 ,nc 0 ,nr ,aa 0))
 261          (merror "amatrix_multiply: inconformable operands")))
 262       (t
 263         `((mnctimes) ,m1 ,m2)))))
 264
 265 (defun dot (n a a0 ainc b b0 binc)
 266   (let ((ai a0) (bi b0) (sum 0))
 267     (dotimes (i n)
 268       (setq
 269         sum (+ sum (* (aref a ai) (aref b bi)))
 270         ai (+ ai ainc)
 271         bi (+ bi binc)))
 272     sum))