src/newdet.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   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 varlist genvar aryp))
  21
  22 ;;these are general type arrays
  23
  24 (defvar *i*)
  25 (defvar *minor1*)
  26 (defvar *binom*)
  27 (defvar *input*)
  28
  29 (defmfun $newdet (mat)
  30   (cond ((not (or (mbagp mat) ($matrixp mat)))
  31          (if ($scalarp mat) mat (list '(%newdet simp) mat)))
  32         (t
  33          (setq mat (check mat))
  34          (unless (= (length mat) (length (cadr mat)))
  35            (merror
  36              (intl:gettext
  37                "newdet: Matrix must be square; found ~M rows, ~M columns.")
  38             (length (cdr mat))
  39             (length (cdadr mat))))
  40          (newdet mat (length (cdr mat)) nil))))
  41
  42 (defmfun $permanent (mat)
  43   (cond ((not (or (mbagp mat) ($matrixp mat)))
  44          (if ($scalarp mat) mat (list '(%permanent simp) mat)))
  45         (t
  46          (setq mat (check mat))
  47          (unless (= (length mat) (length (cadr mat)))
  48            (merror
  49              (intl:gettext
  50                "permanent: Matrix must be square; found ~M rows, ~M columns.")
  51             (length (cdr mat))
  52             (length (cdadr mat))))
  53          (newdet mat (length (cdr mat)) t))))
  54
  55 (defun newdet (a n perm)
  56   (prog (rr k j old new vlist m loc addr sign)
  57      (when (> n 50)
  58        (merror (intl:gettext "newdet: matrix must be 50 by 50 or smaller; found size: ~M") n))
  59      (setq  *binom* (make-array (list (1+ n) (1+ n)) :element-type 'integer))
  60      (setq  *minor1* (make-array (list 2 (1+ (setq rr (pascal n))))))
  61      (setq  *i* (make-array (+ 2 n)))
  62      (do ((k 0 (1+ k)))
  63          ((> k 1))
  64        (do ((j 0 (1+ j)))
  65            ((> j rr))
  66          (setf (aref *minor1* k j) '(0 . 1))))
  67      (do ((k 0 (1+ k)))
  68          ((> k (1+ n)))
  69        (setf (aref *i* k) -1))
  70      (setq  *input* (make-array (list (1+ n) (1+ n))))
  71      (do ((k 1 (1+ k)))
  72          ((> k n))
  73        (do ((j 1 (1+ j)))
  74            ((> j n))
  75          (newvar1 (setf (aref *input* k j) (let ((aryp t)) (maref a k j))))))
  76      (newvar (cons '(mtimes) vlist))
  77      (do ((k 1 (1+ k)))
  78          ((> k n))
  79        (do ((j 1 (1+ j)))
  80            ((> j n))
  81          (setf (aref *input* k j) (cdr (ratrep* (aref *input* k j))))))
  82      (setq new 1)
  83      (setq old 0)
  84      (setf (aref *i* 0) n)
  85      (do ((loc 1 (1+ loc)))
  86          ((> loc n))
  87        (setf (aref *minor1* old (1- loc)) (aref *input* 1 loc)))
  88      (setq m 1)
  89      g0193 (when (> m (1- n)) (go ret))
  90      (setq loc 0)
  91      (setq j 1)
  92      g0189 (when (> j m) (go nextminor))
  93      (setf (aref *i* j) (- m j))
  94      (incf j)
  95      (go g0189)
  96      nextminor
  97      (cond ((not (equal (aref *minor1* old loc) '(0 . 1)))
  98             (setq k (1- n))
  99             (setq j 0)
 100             (setq addr (+ loc (aref *binom* k (1+ m))))
 101             (setq sign 1))
 102            (t (go over)))
 103      nextuse
 104      (cond
 105        ((equal k (aref *i* (1+ j)))
 106         (incf j)
 107         (setq sign (- sign)))
 108        (t
 109         (setf (aref *minor1* new addr)
 110               (ratplus
 111                (aref *minor1* new addr)
 112                (rattimes (aref *minor1* old loc)
 113                          (cond ((or (= sign 1) perm)
 114                                 (aref *input* (1+ m) (1+ k)))
 115                                (t (ratminus (aref *input* (1+ m) (1+ k)))))
 116                          t)))))
 117      (when (> k 0)
 118        (decf k)
 119        (decf addr (aref *binom* k (- m j)))
 120        (go nextuse))
 121      (setf (aref *minor1* old loc)  '(0 . 1))
 122      over (incf loc)
 123      (setq j m)
 124      back (when (> 1 j)
 125             (incf m)
 126             (setq old (- 1 old))
 127             (setq new (- 1 new))
 128             (go g0193))
 129      (setf (aref *i* j) (1+ (aref *i* j)))
 130      (if (> (aref *i* (1- j)) (aref *i* j))
 131          (go nextminor)
 132          (setf (aref *i* j) (- m j)))
 133
 134      (decf j)
 135      (go back)
 136      ret
 137      (return (cons (list 'mrat 'simp varlist genvar) (aref *minor1* old 0)))))
 138
 139 (defun pascal (n)
 140   (setf (aref *binom* 0 0) 1)
 141   (do ((h 1 (1+ h)))
 142       ((> h n) (1- (aref *binom* n (ash n -1))))
 143     (setf (aref *binom* h 0) 1)
 144     (setf (aref *binom* (1- h) h) 0)
 145     (do ((j 1 (1+ j)))
 146         ((> j h))
 147       (setf (aref *binom* h j) (+ (aref *binom* (1- h) (1- j)) (aref *binom* (1- h) j))))))