src/newdet.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module newdet)
  14
  15 ;; THIS IS A VERSION OF THE GENTLEMAN-JOHNSON TREE-MINOR DETERMINANT
  16 ;; USING RATIONAL FUNCTIONS.  "A" CAN BE A MATRIX OR AN ARRAY.
  17 ;; ANSWER IS IN RATIONAL FORM.
  18 ;; RJF  5/2/73
  19
  20 (declare-top (special vlist))
  21
  22 ;;these are general type arrays
  23
  24 (defvar *i*)
  25 (defvar *minor1*)
  26 (defvar *binom*)
  27 (defvar *input*)
  28
  29 #|
  30   Notes by Michel Talon, 2019-02-05.
  31
  32   The algorithm in newdet:
  33   condider computing the determinant of the 4x4 matrix. The essence of the algorithm
  34   is to  store as much intermediary results as possible to avoid recomputation and
  35   group stuff to minimise the number of multiplications. The idea is that when
  36   matrix elements are rationsl functions, multiplication is expensive, and one first
  37   needs to minimize them.
  38
  39                           [ a      a      a      a     ]
  40                           [  1, 1   1, 2   1, 3   1, 4 ]
  41                           [                            ]
  42                           [ a      a      a      a     ]
  43                           [  2, 1   2, 2   2, 3   2, 4 ]
  44   (%o1)                   [                            ]
  45                           [ a      a      a      a     ]
  46                           [  3, 1   3, 2   3, 3   3, 4 ]
  47                           [                            ]
  48                           [ a      a      a      a     ]
  49                           [  4, 1   4, 2   4, 3   4, 4 ]
  50   First one focuses on the first column. One stores the 4 elements
  51   a[1,1],a[2,1],a[3,1],a[4,1]  somewhere.
  52   Then one focuses on the first and second column and one computes the 6 minors
  53   2x2 which are:
  54                                  [ a      a     ]
  55                                  [  1, 1   1, 2 ]
  56   (%o2)                          [              ]
  57                                  [ a      a     ]
  58                                  [  2, 1   2, 2 ]
  59
  60                                  [ a      a     ]
  61                                  [  1, 1   1, 2 ]
  62   (%o3)                          [              ]
  63                                  [ a      a     ]
  64                                  [  3, 1   3, 2 ]
  65   etc. There are 6 =C(4,2) minors because there are 6 ways to choose 2 rows between 4.
  66   Ech minor is computed developing on the last column. So for example the first one is
  67   - a[1,2]*a[2,1] + a[2,2]*a[1,1]
  68   where a[1,1] and a[2,1] are already stored. This takes 2 multiplications hence in total
  69   2*6=12 mukltiplications.
  70   Then one focuses on the columns 1,2,3 and compute and store the 4 minors of type 3x3
  71   such as
  72                               [ a      a      a     ]
  73                               [  1, 1   1, 2   1, 3 ]
  74                               [                     ]
  75   (%o4)                       [ a      a      a     ]
  76                               [  2, 1   2, 2   2, 3 ]
  77                               [                     ]
  78                               [ a      a      a     ]
  79                               [  3, 1   3, 2   3, 3 ]
  80   For each one develop on the last column, which gives 3 multiplications of
  81   a[1,3], a[2,3], a[3,3] by 2x2 minors already computed and stored. In total
  82   3*4=12 multiplications. Finally we develop the original determinant on the last column
  83   giving 4 multiplications by already stored 3x3 minors.
  84
  85   Finally the determinant is computed using 12+12+4=28 multiplications that is
  86   n*(2^(n-1)-1) multiplications with n=4. It is easy to see that this formula is general.
  87   This has to be compared with the straightforward computation of the determinant in which
  88   we have 4! = 24 terms each one requiring 3 multiplications thus 72 multiplications.
  89   The gain comes from the fact that intermediate multiplications are hidden in the minors
  90   and stored by the algorithm. The algorithm has to manage storing the values of the
  91   intermediate minors in arrays of size 4 then 6 then 4 in which one recognizes the fourth row
  92   of the Pascal triangle. It also has to keep track of the alternating + and - signs when
  93   developing on last column. Finally all computations are done using rational algorithms
  94   in maxima in the few following lines:
  95    (rattimes (aref *minor1* old loc)
  96                          (cond ((or (= sign 1) perm)
  97                                 (aref *input* (1+ m) (1+ k)))
  98                                (t (ratminus (aref *input* (1+ m) (1+ k)))))
  99                          t)
 100   where *minor1* keeps the precomputed minors, and rattimes does the product. As to
 101   Gentleman and Johnson paper, its aim is to show the the procedure described above is optimal
 102   with the choice of possible different strategies of grouping.
 103  |#
 104
 105 (defmfun $newdet (mat)
 106   (cond ((not (or (mbagp mat) ($matrixp mat)))
 107          (if ($scalarp mat) mat (list '(%newdet simp) mat)))
 108         (t
 109          (setq mat (check mat))
 110          (unless (= (length mat) (length (cadr mat)))
 111            (merror
 112              (intl:gettext
 113                "newdet: Matrix must be square; found ~M rows, ~M columns.")
 114             (length (cdr mat))
 115             (length (cdadr mat))))
 116          (newdet mat (length (cdr mat)) nil))))
 117
 118 (defmfun $permanent (mat)
 119   (cond ((not (or (mbagp mat) ($matrixp mat)))
 120          (if ($scalarp mat) mat (list '(%permanent simp) mat)))
 121         (t
 122          (setq mat (check mat))
 123          (unless (= (length mat) (length (cadr mat)))
 124            (merror
 125              (intl:gettext
 126                "permanent: Matrix must be square; found ~M rows, ~M columns.")
 127             (length (cdr mat))
 128             (length (cdadr mat))))
 129          (newdet mat (length (cdr mat)) t))))
 130
 131 (defun newdet (a n perm)
 132   (prog (rr k j old new vlist m loc addr sign)
 133      (when (> n 50)
 134        (merror (intl:gettext "newdet: matrix must be 50 by 50 or smaller; found size: ~M") n))
 135      (setq  *binom* (make-array (list (1+ n) (1+ n)) :element-type 'integer))
 136      (setq  *minor1* (make-array (list 2 (1+ (setq rr (pascal n))))))
 137      (setq  *i* (make-array (+ 2 n)))
 138      (do ((k 0 (1+ k)))
 139          ((> k 1))
 140        (do ((j 0 (1+ j)))
 141            ((> j rr))
 142          (setf (aref *minor1* k j) '(0 . 1))))
 143      (do ((k 0 (1+ k)))
 144          ((> k (1+ n)))
 145        (setf (aref *i* k) -1))
 146      (setq  *input* (make-array (list (1+ n) (1+ n))))
 147      (do ((k 1 (1+ k)))
 148          ((> k n))
 149        (do ((j 1 (1+ j)))
 150            ((> j n))
 151          (newvar1 (setf (aref *input* k j) (let ((aryp t)) (maref a k j))))))
 152      (newvar (cons '(mtimes) vlist))
 153      (do ((k 1 (1+ k)))
 154          ((> k n))
 155        (do ((j 1 (1+ j)))
 156            ((> j n))
 157          (setf (aref *input* k j) (cdr (ratrep* (aref *input* k j))))))
 158      (setq new 1)
 159      (setq old 0)
 160      (setf (aref *i* 0) n)
 161      (do ((loc 1 (1+ loc)))
 162          ((> loc n))
 163        (setf (aref *minor1* old (1- loc)) (aref *input* 1 loc)))
 164      (setq m 1)
 165      g0193 (when (> m (1- n)) (go ret))
 166      (setq loc 0)
 167      (setq j 1)
 168      g0189 (when (> j m) (go nextminor))
 169      (setf (aref *i* j) (- m j))
 170      (incf j)
 171      (go g0189)
 172      nextminor
 173      (cond ((not (equal (aref *minor1* old loc) '(0 . 1)))
 174             (setq k (1- n))
 175             (setq j 0)
 176             (setq addr (+ loc (aref *binom* k (1+ m))))
 177             (setq sign 1))
 178            (t (go over)))
 179      nextuse
 180      (cond
 181        ((equal k (aref *i* (1+ j)))
 182         (incf j)
 183         (setq sign (- sign)))
 184        (t
 185         (setf (aref *minor1* new addr)
 186               (ratplus
 187                (aref *minor1* new addr)
 188                (rattimes (aref *minor1* old loc)
 189                          (cond ((or (= sign 1) perm)
 190                                 (aref *input* (1+ m) (1+ k)))
 191                                (t (ratminus (aref *input* (1+ m) (1+ k)))))
 192                          t)))))
 193      (when (> k 0)
 194        (decf k)
 195        (decf addr (aref *binom* k (- m j)))
 196        (go nextuse))
 197      (setf (aref *minor1* old loc)  '(0 . 1))
 198      over (incf loc)
 199      (setq j m)
 200      back (when (> 1 j)
 201             (incf m)
 202             (setq old (- 1 old))
 203             (setq new (- 1 new))
 204             (go g0193))
 205      (setf (aref *i* j) (1+ (aref *i* j)))
 206      (if (> (aref *i* (1- j)) (aref *i* j))
 207          (go nextminor)
 208          (setf (aref *i* j) (- m j)))
 209
 210      (decf j)
 211      (go back)
 212      ret
 213      (return (cons (list 'mrat 'simp varlist genvar) (aref *minor1* old 0)))))
 214
 215 (defun pascal (n)
 216   (setf (aref *binom* 0 0) 1)
 217   (do ((h 1 (1+ h)))
 218       ((> h n) (1- (aref *binom* n (ash n -1))))
 219     (setf (aref *binom* h 0) 1)
 220     (setf (aref *binom* (1- h) h) 0)
 221     (do ((j 1 (1+ j)))
 222         ((> j h))
 223       (setf (aref *binom* h j) (+ (aref *binom* (1- h) (1- j)) (aref *binom* (1- h) j))))))