/usr/share/maxima/5.32.1/src/newdet.lisp is in maxima-src 5.32.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 | ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; The data in this file contains enhancments. ;;;;;
;;; ;;;;;
;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
;;; All rights reserved ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; (c) Copyright 1980 Massachusetts Institute of Technology ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package :maxima)
(macsyma-module newdet)
;; THIS IS A VERSION OF THE GENTLEMAN-JOHNSON TREE-MINOR DETERMINANT
;; USING RATIONAL FUNCTIONS. "A" CAN BE A MATRIX OR AN ARRAY.
;; ANSWER IS IN RATIONAL FORM.
;; RJF 5/2/73
(declare-top (special vlist varlist genvar aryp))
;;these are general type arrays
(defvar *i*)
(defvar *minor1*)
(defvar *binom*)
(defvar *input*)
(defmfun $newdet (mat)
(cond ((not (or (mbagp mat) ($matrixp mat)))
(if ($scalarp mat) mat (list '(%newdet simp) mat)))
(t
(setq mat (check mat))
(unless (= (length mat) (length (cadr mat)))
(merror
(intl:gettext
"newdet: Matrix must be square; found ~M rows, ~M columns.")
(length (cdr mat))
(length (cdadr mat))))
(newdet mat (length (cdr mat)) nil))))
(defmfun $permanent (mat)
(cond ((not (or (mbagp mat) ($matrixp mat)))
(if ($scalarp mat) mat (list '(%permanent simp) mat)))
(t
(setq mat (check mat))
(unless (= (length mat) (length (cadr mat)))
(merror
(intl:gettext
"permanent: Matrix must be square; found ~M rows, ~M columns.")
(length (cdr mat))
(length (cdadr mat))))
(newdet mat (length (cdr mat)) t))))
(defun newdet (a n perm)
(prog (rr k j old new vlist m loc addr sign)
(when (> n 50)
(merror (intl:gettext "newdet: matrix must be 50 by 50 or smaller; found size: ~M") n))
(setq *binom* (make-array (list (1+ n) (1+ n)) :element-type 'integer))
(setq *minor1* (make-array (list 2 (1+ (setq rr (pascal n))))))
(setq *i* (make-array (+ 2 n)))
(do ((k 0 (1+ k)))
((> k 1))
(do ((j 0 (1+ j)))
((> j rr))
(setf (aref *minor1* k j) '(0 . 1))))
(do ((k 0 (1+ k)))
((> k (1+ n)))
(setf (aref *i* k) -1))
(setq *input* (make-array (list (1+ n) (1+ n))))
(do ((k 1 (1+ k)))
((> k n))
(do ((j 1 (1+ j)))
((> j n))
(newvar1 (setf (aref *input* k j) (let ((aryp t)) (maref a k j))))))
(newvar (cons '(mtimes) vlist))
(do ((k 1 (1+ k)))
((> k n))
(do ((j 1 (1+ j)))
((> j n))
(setf (aref *input* k j) (cdr (ratrep* (aref *input* k j))))))
(setq new 1)
(setq old 0)
(setf (aref *i* 0) n)
(do ((loc 1 (1+ loc)))
((> loc n))
(setf (aref *minor1* old (1- loc)) (aref *input* 1 loc)))
(setq m 1)
g0193 (when (> m (1- n)) (go ret))
(setq loc 0)
(setq j 1)
g0189 (when (> j m) (go nextminor))
(setf (aref *i* j) (- m j))
(incf j)
(go g0189)
nextminor
(cond ((not (equal (aref *minor1* old loc) '(0 . 1)))
(setq k (1- n))
(setq j 0)
(setq addr (+ loc (aref *binom* k (1+ m))))
(setq sign 1))
(t (go over)))
nextuse
(cond
((equal k (aref *i* (1+ j)))
(incf j)
(setq sign (- sign)))
(t
(setf (aref *minor1* new addr)
(ratplus
(aref *minor1* new addr)
(rattimes (aref *minor1* old loc)
(cond ((or (= sign 1) perm)
(aref *input* (1+ m) (1+ k)))
(t (ratminus (aref *input* (1+ m) (1+ k)))))
t)))))
(when (> k 0)
(decf k)
(decf addr (aref *binom* k (- m j)))
(go nextuse))
(setf (aref *minor1* old loc) '(0 . 1))
over (incf loc)
(setq j m)
back (when (> 1 j)
(incf m)
(setq old (- 1 old))
(setq new (- 1 new))
(go g0193))
(setf (aref *i* j) (1+ (aref *i* j)))
(if (> (aref *i* (1- j)) (aref *i* j))
(go nextminor)
(setf (aref *i* j) (- m j)))
(decf j)
(go back)
ret
(return (cons (list 'mrat 'simp varlist genvar) (aref *minor1* old 0)))))
(defun pascal (n)
(setf (aref *binom* 0 0) 1)
(do ((h 1 (1+ h)))
((> h n) (1- (aref *binom* n (ash n -1))))
(setf (aref *binom* h 0) 1)
(setf (aref *binom* (1- h) h) 0)
(do ((j 1 (1+ j)))
((> j h))
(setf (aref *binom* h j) (+ (aref *binom* (1- h) (1- j)) (aref *binom* (1- h) j))))))
|