src/math/row-operns.lisp

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