share/vector/vector_rebuild.lisp

   1
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3
   4 ;; Some utilities for working with vectors.
   5
   6 ;; Copyright (C)  Nov. 2008  Volker van Nek
   7
   8 ;; modified 08-12-03: vector_factor factors lists and matrices
   9 ;;          08-12-05: vector_eval: $ratprint set to false
  10 ;;          08-12-10: rename stardisp to stardisp1, assign property of $stardisp
  11 ;;          08-12-14: vector_eval: cut out sratsimp, rename it to vector_simp
  12 ;;                    $vector_rebuild: evaluate mnctimes, bugfix case mdefine
  13 ;;          20-02-27: TeX-support added
  14
  15 ;; This program is free software; you can redistribute it and/or modify
  16 ;; it under the terms of the GNU General Public License as published by
  17 ;; the Free Software Foundation; either version 2 of the License, or
  18 ;; (at your option) any later version.
  19
  20 ;; This program is distributed in the hope that it will be useful,
  21 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  22 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  23 ;; GNU General Public License for more details.
  24
  25 ;; You should have received a copy of the GNU General Public License
  26 ;; along with this program; if not, write to the Free Software
  27 ;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  28 ;; MA 02110-1301, USA.
  29
  30 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  31
  32 ;; Functions at Maxima level:
  33 ;; vector_rebuild(x), vector_rebuild(x,[param_list])
  34 ;; vector_simp(x),    has evfun property
  35 ;; vector_factor(x),  has evfun property
  36 ;; extract_equations(x)
  37
  38 ;; Operators at Maxima level:
  39 ;; |x| vector length
  40 ;; x~y cross product (lbp=134, rbp=133)
  41
  42 ;; Option variable at Maxima level:
  43 ;; vector_factor_minus
  44
  45 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  46
  47 ;; mtimesq
  48
  49 ;; In combination with setting listarith & friends to false
  50 ;; mtimesq allows to suppress automatic arithmetrics
  51 ;; and to bypass the displa function dim-mquotient:
  52 ;;                  1                                       1 [ 2 ]
  53 ;;   1/r*[2,3]; --> - [2, 3]   or   1/r*matrix([2],[3]);--> - [   ]
  54 ;;                  r                                       r [ 3 ]
  55 ;; mtimesq display symbol is "*" and is settable by stardisp.
  56 ;; One evaluation steps from mtimesq to mtimes and vice versa.
  57
  58 (meval `(($infix) mtimesq 120 119))
  59 (defprop mtimesq simp-mtimesq operators)
  60
  61 (defun simp-mtimesq (a tmp z)
  62   (declare (ignore tmp))
  63   `((mtimesq simp) ,(simplifya (cadr a) z) ,(simplifya (caddr a) z)))
  64
  65 (defmfun mtimesq (a b) `((mtimes simp) ,a ,b))
  66
  67 (putprop 'mtimesq (get 'mtimes 'dissym) 'dissym)
  68
  69 ;; extend stardisp in displa.lisp :
  70 (defun stardisp1 (symbol val)
  71   (declare (ignore symbol))
  72   (putprop 'mtimes (if val '(#\*) '(#\space)) 'dissym)
  73   (putprop 'mtimesq (get 'mtimes 'dissym) 'dissym) )
  74 ;;
  75 (defprop $stardisp stardisp1 assign)
  76
  77 (defun mtimesqp (expr)
  78   (and (not (atom expr)) (eq 'mtimesq (caar expr))) )
  79
  80 ;; TeX-support:
  81 (putprop 'mtimesq 'tex-mtimesq 'tex)
  82 ;;
  83 (defun tex-mtimesq (x l r)
  84   (append l
  85           (tex (cadr x) nil nil 'mparen 'mparen)
  86           `("\\,")
  87           ; `("\\cdot") ; if preferred
  88           (tex (caddr x) nil nil 'mparen 'mparen)
  89            r ))
  90
  91 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  92
  93 ;; $vector_simp
  94
  95 (defun $vector_simp (expre$$ion)
  96   (let (($listarith t) ($doallmxops t))
  97     ($expand (meval `(($ev) ,expre$$ion $infeval))) ))
  98     ;; mtimesq needs an extra evaluation here
  99
 100 (putprop '$vector_simp t 'evfun)
 101
 102 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 103
 104 ;; $vector_rebuild
 105
 106 (defun $vector_rebuild (expr &optional (params '((mlist))) )
 107   (when (or (not ($listp params))
 108             (some #'(lambda (p) (or ($constantp p)
 109                                     (not (or ($symbolp p) ($subvarp p)))))
 110                   (cdr params)))
 111     (merror
 112       (format nil "The optional second argument to `vector_rebuild'~%~
 113         must be a list of symbols or subscripted variables.")))
 114   (if (or (atom expr) ($scalarp expr))
 115     expr
 116     (let ((op (mop expr))
 117           (args (margs expr)) )
 118       (if (not (car args)) (return-from $vector_rebuild expr)) ;; noarg-op
 119       (cond
 120         ((eq op 'mequal)
 121           (cons `(mequal simp)
 122             (mapcar #'(lambda (arg) ($vector_rebuild arg params)) args) ))
 123         ((eq op 'mdefine)
 124           (meval `((mdefine simp)
 125             ,(cadr expr)
 126             ,($vector_rebuild (caddr expr) params)) ))
 127         (t
 128           (vector-rebuild expr (cdr params)) )))))
 129
 130 (defun vector-rebuild (expr params)
 131   (let (coef-matrix col-flag tmp-vec col
 132         (res '((mplus) 0))
 133         (n 0)
 134         (vec ($vector_simp expr)) )
 135     (when (not (vector-p vec)) (return-from vector-rebuild expr))
 136     (when (zero-vector-p vec) (return-from vector-rebuild vec))
 137     (when (column-vector-p vec) (setq col-flag t))
 138     (dolist (a (cdr vec))
 139       (when col-flag (setq a (cadr a)))
 140       (setq coef-matrix
 141         (append coef-matrix (list (coef-list a params)))) )
 142     (setq params (append params '(1)))
 143     (dolist (p params)
 144       (setq col (mapcar #'(lambda (row) (nth n row)) coef-matrix))
 145       (incf n)
 146       (when (some #'(lambda (e) (not (zerop1 e))) col)
 147         (progn
 148           (setq tmp-vec `((mlist simp) ,@col) )
 149           (when col-flag
 150             (setq tmp-vec (row-to-column tmp-vec)))
 151           (setq res
 152             (append res `(((mtimes) ,p ,tmp-vec)) ))) ))
 153     (meval res) ))
 154
 155 ;; e.g. expr=4-2*t+3*s, params=(s t) --> (3 -2 4)
 156 (defun coef-list (expr params)
 157   (let (res c)
 158     (dolist (p params)
 159       (setq c (meval `(($coeff) ,expr ,p)))
 160       (push c res)
 161       (setq expr
 162         (meval `((mplus) ,expr ((mtimes) ((mminus) ,c) ,p))) ))
 163     (reverse (cons expr res)) ))
 164
 165 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 166
 167 ;; $extract_equations
 168
 169 (defun $extract_equations (equation)
 170   (let ((err-msg "Argument to `extract_equations' is not a vector equation.")
 171         args left right)
 172     (and (not (atom equation))
 173          (not (eq (mop equation) 'mequal))
 174          (merror err-msg))
 175     (setq args (cdr equation))
 176     (setq left ($vector_simp (car args)))
 177     (setq right ($vector_simp (cadr args)))
 178     (cond
 179       ((and (vector-p left) (vector-p right))) ;; OK
 180       ;; due to a bug in Maxima allow left or right to be zero:
 181       ((and (eql 0 left) (vector-p right))
 182         (setq left
 183           (cons '(mlist simp)
 184             (make-list (vector-dim right) :initial-element 0))) )
 185       ((and (eql 0 right) (vector-p left))
 186         (setq right
 187           (cons '(mlist simp)
 188             (make-list (vector-dim left) :initial-element 0))) )
 189       (t
 190          (merror err-msg) ))
 191     (when (not ($listp left)) (setq left (column-to-row left)))
 192     (when (not ($listp right)) (setq right (column-to-row right)))
 193     (meval `(($map) "=" ,left ,right)) ))
 194
 195 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 196
 197 ;; $vector_factor
 198
 199 (defun $vector_factor (expr)
 200   (if (atom expr)
 201     expr
 202     (let ((op (mop expr)))
 203       (cond
 204         ((eq op 'mequal)
 205           `((mequal simp) ,@(mapcar #'$vector_factor (cdr expr))) )
 206         ((eq op 'mdefine)
 207           (meval `((mdefine simp)
 208             ,(cadr expr)
 209             ,($vector_factor (caddr expr))) ))
 210         ((member op '(mplus mnctimes crossq) :test #'eq)
 211           `((,op simp) ,@(mapcar #'$vector_factor (cdr expr))) )
 212         (t
 213           (vector-extract-gcd expr) ))))) ;; incl. op=mtimes,mminus
 214
 215 (putprop '$vector_factor t 'evfun)
 216
 217 ;; if you wish to factor out a minus sign, set this to true
 218 (defmvar $vector_factor_minus nil)
 219
 220 (defun vector-extract-gcd (expr)
 221   (let (fac vec args minus-flag $ratprint)
 222     (setq vec ($vector_simp expr))
 223     (cond
 224       ((or (zero-vector-p vec) (zero-$matrix-p vec))
 225         vec)
 226       ((or ($listp vec) ($matrixp vec))
 227         (setq args
 228           (if ($listp vec)
 229             (cdr vec)
 230             (apply #'append (mapcar #'cdr (cdr vec))) ))
 231         (setq fac
 232           (reduce #'(lambda (a b) ($gcd a b)) (cons (car args) args)))
 233         (setq vec
 234           (meval
 235             `(($fullmapl)
 236               ((lambda) ((mlist) e) (($ratsimp) ((mquotient) e ,fac)))
 237               ,vec )))
 238         (and
 239           $vector_factor_minus
 240           (every
 241             #'(lambda (a) (eq t (meval `(($is) ((mlessp) ,a 0)))))
 242             args)
 243           (setq minus-flag t)
 244           (setq vec
 245             (meval
 246               `(($fullmapl) ((lambda) ((mlist) e) ((mminus) e)) ,vec) )))
 247         (if minus-flag
 248           (if (eql 1 fac)
 249             `((mminus) ,vec)
 250             `((mminus) ((mtimesq) ,fac ,vec)) )
 251           (if (eql 1 fac)
 252             vec
 253             `((mtimesq) ,fac ,vec) )))
 254       (t
 255         expr ))))
 256
 257 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 258
 259 ;; vector length
 260
 261 (meval `(($matchfix) $\| $\|))
 262
 263 (defprop $\| simp-vector-length operators)
 264
 265 (defun simp-vector-length (a tmp z)
 266   (declare (ignore tmp) (ignore z))
 267   (oneargcheck a)
 268   (setq a ($vector_simp (cadr a)))
 269   (cond
 270     (($scalarp a)
 271       (meval `((mabs simp) ,a)) )
 272     ((vector-p a)
 273       (let (args)
 274         (setq args
 275           (if ($listp a)
 276             (cdr a)
 277             (mapcar #'cadr (cdr a)) ))
 278         (setq args (mapcar #'(lambda (a) (meval `((mexpt) ,a 2))) args))
 279         (meval `(($sqrt simp) ((mplus simp) ,@args))) ))
 280     (t
 281       `(($\|) ,a) )))
 282
 283 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 284
 285 ;; cross product
 286
 287 (meval `(($infix) $\~ 134 133)) ;; lbp and rbp copied from vect.mac
 288
 289 (defprop $\~ simp-vector-cross operators)
 290
 291 ;; if it seems fit we can suppress automatic simplification
 292 ;; by using an alias:
 293 ;;
 294 ;; (meval `(($infix) crossq 134 133))
 295 ;; (putprop 'crossq '(#\~) 'dissym)
 296 ;;
 297 ;; (defun simp-vector-cross (a tmp z)
 298 ;;   (declare (ignore tmp) (ignore z))
 299 ;;   (twoargcheck a)
 300 ;;   `((crossq simp) ,(cadr a) ,(caddr a)) )
 301 ;;    ;; again, here is one more step of evaluation needed
 302
 303 ;; (defmfun crossq (a b)
 304 ;;   (let ((v ($vector_simp a))
 305 ;;         (w ($vector_simp b))
 306 ;;         cross col-flag dim-v dim-w) ;;))
 307
 308 (defun simp-vector-cross (a tmp z)
 309   (declare (ignore tmp) (ignore z))
 310   (twoargcheck a)
 311   (let ((v ($vector_simp (cadr a)))
 312         (w ($vector_simp (caddr a)))
 313         cross col-flag dim-v dim-w)
 314     (if (and (vector-p v) (vector-p w))
 315       (progn
 316         (when (column-vector-p v)
 317           (setq v (column-to-row v) col-flag t))
 318         (when (column-vector-p w)
 319           (setq w (column-to-row w) col-flag t))
 320         (setq dim-v (1- (length v)) dim-w (1- (length w)))
 321         (cond
 322           ((and (= 2 dim-v) (= 2 dim-w))
 323             ;; v[1]*w[2]-v[2]*w[1]
 324             (meval `((mplus simp)
 325               ((mtimes simp) ,(cadr v) ,(caddr w))
 326               ((mminus) ((mtimes simp) ,(caddr v) ,(cadr w))) )))
 327           ((and (= 3 dim-v) (= 3 dim-w))
 328             ;; [ v[2]*w[3]-v[3]*w[2],
 329             ;;   v[3]*w[1]-v[1]*w[3],
 330             ;;   v[1]*w[2]-v[2]*w[1] ]
 331             (setq cross
 332               `((mlist simp)
 333                 ,(meval `((mplus simp)
 334                   ((mtimes simp) ,(caddr v) ,(cadddr w))
 335                   ((mminus) ((mtimes simp) ,(cadddr v) ,(caddr w))) ))
 336                 ,(meval `((mplus simp)
 337                   ((mtimes simp) ,(cadddr v) ,(cadr w))
 338                   ((mminus) ((mtimes simp) ,(cadr v) ,(cadddr w))) ))
 339                 ,(meval `((mplus simp)
 340                   ((mtimes simp) ,(cadr v) ,(caddr w))
 341                   ((mminus) ((mtimes simp) ,(caddr v) ,(cadr w))) )) ))
 342             (if col-flag
 343               (row-to-column cross)
 344               cross ))
 345           (t
 346             `(($\~) ,v ,w) ) ))
 347       `(($\~) ,v ,w) )))
 348
 349 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 350
 351 ;; some small helpers
 352
 353 (defun column-to-row (col) ;; assume col to be a column vector (a $matrix)
 354   (cons '(mlist simp)
 355     (mapcar #'cadr (cdr col)) ))
 356
 357 (defun row-to-column (row) ;; assume row to be a row vector (a mlist)
 358   (cons '($matrix simp)
 359     (mapcar #'(lambda (e) (list '(mlist simp) e))
 360       (cdr row) )))
 361
 362 (defun vector-dim (vec) ;; assume vec to be a row or column vector
 363   (length (cdr vec)) )
 364
 365 (defun column-vector-p (obj)
 366   (and ($matrixp obj)
 367        (= 2 (length (cadr obj))) ))
 368
 369 (defun vector-p (obj)
 370   (or ($listp obj) (column-vector-p obj)))
 371
 372 (defun zero-vector-p (obj)
 373   (or (zero-mlist-p obj)
 374       (and (zero-$matrix-p obj) (= 2 (length (cadr obj)))) ))
 375
 376 (defun zero-mlist-p (obj)
 377   (and ($listp obj) (every #'zerop1 (cdr obj)) ))
 378
 379 (defun zero-$matrix-p (obj)
 380   (and ($matrixp obj) (every #'zero-mlist-p (cdr obj)) ))
 381
 382 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;