src/homog.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 homog)
  14
  15 (load-macsyma-macros ratmac)
  16
  17 (declare-top (special *hvar *hmat))
  18
  19 (defun addvardeg (n l lt)
  20   (mapc #'(lambda (x) (push (cons n x) lt)) l) lt)
  21
  22 (defun ptermvec (p)
  23   (ltermvec p (sort (listovars p) #'pointergp) nil))
  24
  25 (defun nzeros (n l)
  26   (dotimes (i n l)
  27     (push 0 l)))
  28
  29 (defun ltermvec (p vl coef?)
  30   (cond ((null vl) (list (if coef? p nil)))
  31         ((pcoefp p) (list (nzeros (length vl) (if coef? p nil))))
  32         ((pointergp (car vl) (car p))
  33          (addvardeg 0 (ltermvec p (cdr vl) coef?) nil))
  34         (t (do ((p (cdr p) (cddr p))
  35                 (lt nil (addvardeg (car p) (ltermvec (cadr p) (cdr vl)
  36                                                      coef?)
  37                                    lt)))
  38                ((null p) lt)))))
  39
  40 ;car(lv) = list of dependent equations
  41 ;caddr (lv) = correspondence between new columns and old ones.
  42
  43 (defun hlinsolve (mat)
  44   (let ((n (1- (length mat))) (m (length (car mat))) arr ndepvar
  45         (mat (mapcar #'(lambda (x) (mapcar #'- (car mat) x)) (cdr mat))))
  46     (setq arr (make-array (list (1+ (max m n)) (+ 2 m))))
  47     (do ((ml mat (cdr ml))              ;solving for m vars
  48          (i 1 (1+ i)))
  49         ((null ml))
  50       (do ((l (car ml) (cdr l))
  51            (j 1 (1+ j)))
  52           ((null l) (setf (aref arr i j) 0))
  53         (setf (aref arr i j) (car l))))
  54     (setq mat (tfgeli1 arr n (1+ m)))
  55     (and (cadr mat) (merror "HLINSOLVE: inconsistent equations.")) ;shouldn't happen
  56                                         ; # indep equations = n - (car mat)
  57                                         ; # dependent vars = # indep equations
  58                                         ; # indep vars = m - # dependent vars
  59     (setq ndepvar (- n (length (car mat))))
  60     (do ((i (1+ n) (1+ i))) ((> i m))
  61       (do ((j (1+ ndepvar) (1+ j))) ((> j m))
  62         (setf (aref arr i j) 0)))
  63     (do ((i (1+ ndepvar) (1+ i))
  64          (det (abs (aref arr 1 1))))
  65         ((> i m))
  66       (setf (aref arr i i) det))
  67     (cond ((signp g (aref arr 1 1))
  68            (do ((i 1 (1+ i))) ((> i ndepvar))
  69              (do ((j (1+ ndepvar) (1+ j))) ((> j m))
  70                (setf (aref arr i j) (- (aref arr i j)))))))
  71     (do ((l (caddr mat) (cdr l))        ;invert var permutation
  72          (i 1 (1+ i)))
  73         ((null l))
  74       (setf (aref arr 0 (car l)) i))
  75     (do ((varord (caddr mat) (cdr varord))
  76          (i 1 (1+ i)))
  77         ((> i ndepvar)
  78          (do ((varord varord (cdr varord))
  79               (i i (1+ i)))
  80              ((> i m)
  81               (do ((ans)
  82                    (i m (1- i)))
  83                   ((< i 1) ans)
  84                 (push (aref arr i 0) ans)))
  85            (do ((vecl)
  86                 (j m (1- j)))
  87                ((< j 1)
  88                 (setf (aref arr (car varord) 0) vecl))
  89              (push (aref arr (aref arr 0 j) i) vecl))))
  90       (do ((gcd 0)
  91            (j i (1+ j)))
  92           ((or (= gcd 1) (> j m))
  93            (setf (aref arr (car varord) 0)
  94                  (abs (truncate (aref arr i i) gcd))))
  95         (setq gcd (gcd gcd (aref arr i j)))))))
  96 ; returns (mixed list of <reduced exp> and <basis vector for null space>)
  97 ; <reduced expon> corresponds to dependent var
  98 ; <basis vector> corresponds to independent var
  99
 100 (defun hreduce (p &optional (vl (setq *hvar (sort (listovars p) 'pointergp)))
 101                       (hl (setq *hmat (hlinsolve (ltermvec p *hvar nil)))))
 102        (cond ((pcoefp p) p)
 103              ((pointergp (car vl) (car p))
 104               (hreduce p (cdr vl) (cdr hl)))
 105              ((numberp (car hl))
 106               (cons (car p)
 107                     (do ((p (cdr p) (cddr p))
 108                          (red (car hl))
 109                          (ans))
 110                         ((null p) (nreverse ans))
 111                         (push (truncate (car p) red) ans)
 112                         (push (hreduce (cadr p) (cdr vl) (cdr hl)) ans))))
 113              (t (do ((p (cdddr p) (cddr p))
 114                      (sum (hreduce (caddr p) (cdr vl) (cdr hl))
 115                           (pplus sum (hreduce (cadr p) (cdr vl) (cdr hl)))))
 116                     ((null p) sum)))))
 117
 118 (defun hexpand (p &optional (hl *hmat) (vl *hvar))
 119   (if (every #'onep hl)
 120       p
 121       (progn
 122         (do ((hl hl (cdr hl))
 123              (i 1 (1+ i))
 124              (pl (ltermvec p vl t)))
 125             ((null hl) (setq p pl))
 126           (if (and (numberp (car hl)) (not (onep (car hl))))
 127               (do ((pl pl (cdr pl))) ((null pl))
 128                 (do ((term (car pl) (cdr term))
 129                      (j (1- i) (1- j)))
 130                     ((= j 0) (rplaca term (* (car term) (car hl))))))))
 131         (do ((hl hl (cdr hl))
 132              (i 1 (1+ i))
 133              (wtlist)
 134              (newwt))
 135             ((null hl) (hsimp p vl))
 136           (unless (numberp (car hl))
 137             (setq wtlist (mapcar #'(lambda (x) (hdot (car hl) x)) p))
 138             (setq newwt (nth (1- i) (car hl)))
 139             (do ((maxwt (apply #'max wtlist))
 140                  (pl p (cdr pl))
 141                  (wtlist wtlist (cdr wtlist)))
 142                 ((null pl))
 143               (do ((term (car pl) (cdr term))
 144                    (j (1- i) (1- j)))
 145                   ((= j 0) (rplaca term (truncate (- maxwt (car wtlist)) newwt))))))))))
 146
 147 (defun hdot (ht pt)
 148        (do ((ht (cdr ht) (cdr ht))
 149             (pt (cdr pt) (cdr pt))
 150             (sum (* (car ht) (car pt)) (+ sum (* (car ht) (car pt)))))
 151            ((null ht) sum)))
 152
 153 (defun hsimp (pl vl)
 154        (do ((pl (cdr pl) (cdr pl))
 155             (p (hsimp1 (car pl) vl) (pplus p (hsimp1 (car pl) vl))))
 156            ((null pl) p)))
 157
 158 (defun hsimp1 (tl vl)
 159        (cond ((null vl) tl)
 160              ((> (car tl) 0)
 161               (list (car vl) (car tl) (hsimp1 (cdr tl) (cdr vl))))
 162              (t (hsimp1 (cdr tl) (cdr vl)))))