home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / homog.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  5.5 KB  |  165 lines
  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. (macsyma-module homog)
  13.  
  14. (LOAD-MACSYMA-MACROS RATMAC)
  15.  
  16. (declare-top (special *hvar *hmat)
  17.      (*lexpr hreduce hexpand))
  18.  
  19. (defun addvardeg (n l lt)
  20.        (mapc (fn (x) (push (cons n x) lt)) l)
  21.        lt)
  22.  
  23. (defun ptermvec (p)
  24.        (ltermvec p (sort (listovars p) #'pointergp) nil))
  25.  
  26. (defun ltermvec (p vl coef?)
  27.        (cond ((null vl) (list (if coef? p nil)))
  28.          ((pcoefp p) (list (nzeros (length vl) (if coef? p nil))))
  29.          ((pointergp (car vl) (car p))
  30.           (addvardeg 0 (ltermvec p (cdr vl) coef?) nil))
  31.          (t (do ((p (cdr p) (cddr p))
  32.              (lt nil (addvardeg (car p) (ltermvec (cadr p) (cdr vl)
  33.                               coef?)
  34.                     lt)))
  35.             ((null p) lt)))))
  36.  
  37. ;car(lv) = list of dependent equations
  38. ;caddr (lv) = correspondence between new columns and old ones.
  39.  
  40. (defun hlinsolve (mat)
  41.        (let ((n (f1- (length mat))) (m (length (car mat))) arr ndepvar
  42.           (mat (mapcar (fn (x) (mapcar '- (car mat) x)) (cdr mat))))
  43.          (setq arr (array nil t (f1+ (max m n)) (f+ 2 m)))
  44.          (do ((ml mat (cdr ml))            ;solving for m vars
  45.           (i 1 (f1+ i)))
  46.          ((null ml))
  47.          (do ((l (car ml) (cdr l))
  48.               (j 1 (f1+ j)))
  49.              ((null l) (store (arraycall t arr i j) 0))
  50.              (store (arraycall t arr i j) (car l))))
  51.          (setq mat (tfgeli1 arr n (f1+ m)))
  52.          (and (cadr mat) (merror "Inconsistent"))    ;shouldn't happen
  53. ; # indep equations = n - (car mat)
  54. ; # dependent vars = # indep equations
  55. ; # indep vars = m - # dependent vars
  56.          (setq ndepvar (f- n (length (car mat))))
  57.          (do ((i (f1+ n) (f1+ i))) ((> i m))
  58.          (do ((j (f1+ ndepvar) (f1+ j))) ((> j m))
  59.              (store (arraycall t arr i j) 0)))
  60.          (do ((i (f1+ ndepvar) (f1+ i))
  61.           (det (abs (arraycall t arr 1 1))))
  62.          ((> i m))
  63.          (store (arraycall t arr i i) det))
  64.          (cond ((signp g (arraycall t arr 1 1))
  65.             (do ((i 1 (f1+ i))) ((> i ndepvar))
  66.             (do ((j (f1+ ndepvar) (f1+ j))) ((> j m))
  67.                 (store (arraycall t arr i j)
  68.                    (f- (arraycall t arr i j)))))))
  69.          (do ((l (caddr mat) (cdr l))        ;invert var permutation
  70.           (i 1 (f1+ i)))
  71.          ((null l))
  72.          (store (arraycall t arr 0 (car l)) i))
  73.          (do ((varord (caddr mat) (cdr varord))
  74.           (i 1 (f1+ i)))
  75.          ((> i ndepvar)
  76.           (do ((varord varord (cdr varord))
  77.                (i i (f1+ i)))
  78.               ((> i m)
  79.                (do ((ans)
  80.                 (i m (f1- i)))
  81.                ((< i 1) (*rearray arr) ans)
  82.                (push (arraycall t arr i 0) ans)))
  83.               (do ((vecl)
  84.                (j m (f1- j)))
  85.               ((< j 1)
  86.                (store (arraycall t arr (car varord) 0) vecl))
  87.               (push (arraycall t arr
  88.                        (arraycall t arr 0 j) i) vecl))))
  89.          (do ((gcd 0)
  90.               (j i (f1+ j)))
  91.              ((or (= gcd 1) (> j m))
  92.               (store (arraycall t arr (car varord) 0)
  93.                  (abs (// (arraycall t arr i i) gcd))))
  94.              (setq gcd (gcd gcd (arraycall t arr i j)))))
  95.          ))                        
  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 (// (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 (andmapc #'onep hl) p
  120.     (progn
  121.        (do ((hl hl (cdr hl))
  122.        (i 1 (f1+ i))
  123.        (pl (ltermvec p vl t)))
  124.       ((null hl) (setq p pl))
  125.     (if (and (numberp (car hl)) (not (onep (car hl))))
  126.         (do ((pl pl (cdr pl))) ((null pl))
  127.         (do ((term (car pl) (cdr term))
  128.              (j (f1- i) (f1- j)))
  129.             ((= j 0) (rplaca term (f* (car term) (car hl))))))))
  130.       (do ((hl hl (cdr hl))
  131.        (i 1 (f1+ i))
  132.        (wtlist)
  133.        (newwt))
  134.       ((null hl) (hsimp p vl))
  135.     (unless (numberp (car hl))
  136.         (setq wtlist (mapcar (fn (x) (hdot (car hl) x)) p))
  137.         (setq newwt (nth (f1- i) (car hl)))
  138.         (do ((maxwt (apply 'max wtlist))
  139.              (pl p (cdr pl))
  140.              (wtlist wtlist (cdr wtlist)))
  141.             ((null pl))
  142.           (do ((term (car pl) (cdr term))
  143.                (j (f1- i) (f1- j)))
  144.               ((= j 0) (rplaca term (// (f- maxwt (car wtlist))
  145.                         newwt))))))))))
  146.  
  147. (defun hdot (ht pt)
  148.        (do ((ht (cdr ht) (cdr ht))
  149.         (pt (cdr pt) (cdr pt))
  150.         (sum (f* (car ht) (car pt))
  151.          (f+ sum (f* (car ht) (car pt)))))
  152.        ((null ht) sum)))
  153.  
  154. (defun hsimp (pl vl)
  155.        (do ((pl (cdr pl) (cdr pl))
  156.         (p (hsimp1 (car pl) vl)
  157.            (pplus p (hsimp1 (car pl) vl))))
  158.        ((null pl) p)))
  159.  
  160. (defun hsimp1 (tl vl)
  161.        (cond ((null vl) tl)
  162.          ((> (car tl) 0)
  163.           (list (car vl) (car tl) (hsimp1 (cdr tl) (cdr vl))))
  164.          (t (hsimp1 (cdr tl) (cdr vl)))))
  165.