src/math/row-operns.lisp

   1
   2 (in-package :lineal.math)
   3
   4 ; u - kv
   5 (defun tuple-multiple-nsubtrn (u k v)
   6   "Vector /u/ minus scalar /k/ times vector /v/"
   7   (do ((u-rem u (cdr u-rem))
   8        (v-rem v (cdr v-rem)))
   9     ((not (or u-rem v-rem)) u)
  10     (decf (car u-rem) (* k (car v-rem)))))
  11
  12 (defun row-clear-leading-with (u rows)
  13   (declare (inline tuple-multiple-nsubtrn))
  14   "Clear the first column with the leading-1 row."
  15   (loop :for to-clear :in rows :collect
  16         (cdr (tuple-multiple-nsubtrn
  17                to-clear (car to-clear) u))))
  18
  19 (defun row-reduction (a incheck-fn nolead-fn)
  20   (when (and a (car a))
  21     (if (zerop (caar a))
  22       (loop
  23         :with toprow = (car a)
  24         :for rows :on a
  25         :unless (zerop (caar rows)) :do
  26         (progn
  27           (rotatef (caar rows) (car toprow))
  28           (rotatef (cdar rows) (cdr toprow))
  29           (funcall incheck-fn (cdr nat-zero)
  30                    toprow rows)
  31           (return))
  32         :collect (cdar rows) :into nat-zero
  33         :finally (funcall nolead-fn nat-zero))
  34       (funcall incheck-fn nil (car a) (cdr a)))))
  35
  36 (defun det (a &optional (destroy nil)
  37               &aux (curr-det 1))
  38   "Determinant"
  39   (labels
  40     ((fix-in-lead
  41        (row toclear)
  42        (setq curr-det (* curr-det (car row)))
  43        (row-clear-leading-with
  44          (tuple-scalar-ndivisn row (car row))
  45          toclear))
  46      (det-recurse
  47        (tmp-a)
  48        (row-reduction
  49          (setq a tmp-a) #'incheck-fn #'nolead-fn))
  50      (incheck-fn
  51        (nat-zero lead toclear)
  52        (when nat-zero (setq curr-det (- curr-det)))
  53        (det-recurse
  54          (nconc nat-zero
  55                 (fix-in-lead lead toclear))))
  56      (nolead-fn
  57        (nat-zero)
  58        (declare (ignore nat-zero))
  59        (setq curr-det 0)))
  60     (det-recurse (if destroy a (copy-tree a)))
  61     ;^ First call the recursion...^
  62     (if a 0 curr-det)));->... then return the result.
  63
  64 (defun r-ef
  65   (a &optional (destroy nil)
  66      &aux (orig-mtrix (if destroy a (copy-tree a))))
  67   "Row-echelon form"
  68   (labels
  69     ((r-ef-recurse
  70        (tmp-a)
  71        (row-reduction
  72          (setq a tmp-a) #'incheck-fn #'nolead-fn))
  73      (incheck-fn
  74        (nat-zero lead toclear)
  75        (r-ef-recurse
  76          (nconc
  77            nat-zero
  78            (row-clear-leading-with
  79              (tuple-scalar-ndivisn lead (car lead))
  80              toclear))))
  81      (nolead-fn (nat-zero) (r-ef-recurse nat-zero)))
  82     (r-ef-recurse orig-mtrix)
  83     orig-mtrix))
  84
  85 (defun rr-ef
  86   (a &optional (destroy nil)
  87      &aux (also-clear nil)
  88      (orig-mtrix (if destroy a (copy-tree a))))
  89   "Reduced row-echelon form"
  90   (labels
  91     ((rr-ef-recurse
  92        (tmp-a)
  93        (row-reduction
  94          (setq a tmp-a) #'incheck-fn #'nolead-fn))
  95      (incheck-fn
  96        (nat-zero lead toclear)
  97        (tuple-scalar-ndivisn lead (car lead))
  98        (setq also-clear
  99              (cons (cdr lead)
 100                    (row-clear-leading-with
 101                      lead also-clear)))
 102        (rr-ef-recurse
 103          (nconc nat-zero
 104                 (row-clear-leading-with
 105                   lead toclear))))
 106      (nolead-fn
 107        (nat-zero)
 108        ;V No leading 1 can be obtained, go to V
 109        ;V next column in /also-clear/ without V
 110        ;V clearing the current one.           V
 111        (do ((toclear also-clear (cdr toclear)))
 112          ((endp toclear))
 113          (rplaca toclear (cdar toclear)))
 114        (rr-ef-recurse nat-zero)))
 115     (rr-ef-recurse orig-mtrix)
 116     orig-mtrix))
 117
 118 (defun nullspace
 119   (a &key (destroy nil)
 120      ((dimdom dimdom-left) (length (car a)))
 121      &aux (also-clear nil)
 122      ;V Collect basis for nullspace.V
 123      (null-basis nil) (nullity-sofar 0))
 124   "A basis for the nullspace of a matrix."
 125   (labels
 126     ((add-to-basis
 127        ()
 128        (do ((toclear also-clear (cdr toclear))
 129             (newcol
 130               (nconc
 131                 (make-list nullity-sofar
 132                            :initial-element 0)
 133                 (cons -1 (make-list (decf dimdom-left)
 134                                     :initial-element 0)))))
 135          ((endp toclear)
 136           (incf nullity-sofar)
 137           (push newcol null-basis))
 138          (push (pop (car toclear)) newcol)))
 139      (nullspace-recurse
 140        (tmp-a)
 141        (row-reduction
 142          (setq a tmp-a) #'incheck-fn #'nolead-fn))
 143      (incheck-fn
 144        (nat-zero lead toclear)
 145        (tuple-scalar-ndivisn lead (car lead))
 146        (setq also-clear
 147              (cons (cdr lead)
 148                    (row-clear-leading-with
 149                      lead also-clear)))
 150        (decf dimdom-left)
 151        (nullspace-recurse
 152          (nconc nat-zero
 153                 (row-clear-leading-with
 154                   lead toclear))))
 155      (nolead-fn
 156        (nat-zero)
 157        (add-to-basis)
 158        (nullspace-recurse nat-zero)))
 159     (nullspace-recurse (if destroy a (copy-tree a)))
 160     ;V When done with rr-ef, add all  V
 161     ;V remaining columns to the basis.V
 162     (dotimes (iters dimdom-left) (add-to-basis))
 163     null-basis))
 164
 165 (defun mtrix-square-identity (n)
 166   "Create a square identity matrix with /n/ rows and columns."
 167   (loop
 168     :with x = 0
 169     :until (= x n) :collect
 170     (nconc (make-list x :initial-element 0)
 171            (cons 1 (make-list
 172                      (- n (incf x))
 173                      :initial-element 0)))))
 174
 175 (defun mtrix-augment (a b)
 176   "Augment matrix /a/ with /b/"
 177   (let ((tails
 178           (mapcar
 179             #'(lambda (arow brow)
 180                 (rplacd (last arow) brow))
 181             a b)))
 182     ;V Return an "undo" procedure.V
 183     #'(lambda ()
 184         (loop :for atail :in tails
 185               :do (rplacd atail nil)))))
 186
 187 ; todo: make row swaps more efficient,
 188 ; perhaps a full-fledged row-reduction function is needed
 189 (defun mtrix-inverse
 190   (a &optional (row-count (length a)))
 191   "Inverse of a given matrix."
 192   (setq a (copy-tree a))
 193   (mapc #'nconc a (mtrix-square-identity row-count))
 194   (mapcar #'(lambda (row) (nthcdr row-count row))
 195           (rr-ef a)))
 196
 197