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 / csimp2.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  19.9 KB  |  614 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 1982 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module csimp2)
  13.  
  14. (load-macsyma-macros rzmac)
  15.  
  16. (declare-top (GENPREFIX C/#))
  17.  
  18. (declare-top (SPLITFILE PLOG)
  19.      (SPECIAL VAR %P%I VARLIST PLOGABS HALF%PI NN* DN*))
  20.  
  21. (DEFMFUN SIMPPLOG (X VESTIGIAL Z)
  22.        VESTIGIAL ;Ignored.
  23.        (PROG (VARLIST DD CHECK Y)
  24.          (ONEARGCHECK X)
  25.          (SETQ CHECK X)
  26.          (SETQ X (SIMPCHECK (CADR X) Z))
  27.          (COND ((EQUAL 0 X) (MERROR "PLOG(0) is undefined"))
  28.            ((AMONG VAR X) ;This is used in DEFINT. 1/19/81. -JIM
  29.             (RETURN (EQTEST (LIST '(%PLOG) X) CHECK))))
  30.          (NEWVAR X)
  31.          (COND
  32.           ((AND (EQ '$%I (CAR VARLIST)) (NOT (ORMAPC #'ATOM (CDR VARLIST))))
  33.            (SETQ DD (TRISPLIT X))
  34.            (COND ((SETQ Z (PATAN (CAR DD) (CDR DD)))
  35.               (RETURN (ADD2* (SIMPLN (LIST '(%LOG) 
  36.     (SIMPEXPT (LIST '(MEXPT) ($EXPAND (LIST '(MPLUS)
  37.         (LIST '(MEXPT) (CAR DD) 2)
  38.         (LIST '(MEXPT) (CDR DD) 2))) '((RAT) 1 2)) 1 NIL)) 1 T)
  39.                      (LIST '(MTIMES) Z '$%I))))))
  40.           ((AND (FREE X '$%I) (EQ ($SIGN X) '$PNZ))
  41.            (RETURN (EQTEST (LIST '(%PLOG) X) CHECK)))
  42.           ((AND (EQUAL ($IMAGPART X) 0) (SETQ Y ($ASKSIGN X)))
  43.            (COND ((EQ Y '$POS) (RETURN (SIMPLN (LIST '(%LOG) X) 1 T)))
  44.              ((AND PLOGABS (EQ Y '$NEG))
  45.               (RETURN (SIMPLN (LIST '(%LOG) (LIST '(MTIMES) -1 X)) 1 NIL)))
  46.               ((EQ Y '$NEG)
  47.                (RETURN (ADD2 %P%I
  48.                      (SIMPLN (LIST '(%LOG) (LIST '(MTIMES) -1 X)) 1 NIL))))
  49.              (T (MERROR "PLOG(0) is undefined"))))
  50.           ((AND (EQUAL ($IMAGPART (SETQ Z (DIV* X '$%I))) 0)
  51.             (SETQ Y ($ASKSIGN Z)))
  52.            (COND
  53.         ((EQUAL Y '$ZERO) (MERROR "PLOG(0) is undefined"))
  54.         (T (COND ((EQ Y '$POS) (SETQ Y 1))
  55.              ((EQ Y '$NEG) (SETQ Y -1)))
  56.            (RETURN (ADD2* (SIMPLN (LIST '(%LOG)
  57.                         (LIST '(MTIMES) Y Z)) 1 NIL)
  58.                   (LIST '(MTIMES) Y '((RAT) 1 2) '$%I '$%PI)))))))
  59.          (RETURN (EQTEST (LIST '(%PLOG) X) CHECK))))
  60.  
  61. (DEFUN PATAN (R I)
  62.   (LET (($NUMER $NUMER))
  63.        (PROG (A B VAR) 
  64.          (SETQ I (SIMPLIFYA I NIL) R (SIMPLIFYA R NIL))
  65.          (COND ((ZEROP1 R)
  66.             (IF (FLOATP I) (SETQ $NUMER T))
  67.             (SETQ I ($ASKSIGN I))
  68.             (COND ((EQUAL I '$POS) (RETURN (SIMPLIFY HALF%PI)))
  69.               ((EQUAL I '$NEG)
  70.                (RETURN (MUL2 -1 (SIMPLIFY HALF%PI))))
  71.               (T (MERROR "ATAN(0//0) has been generated."))))
  72.            ((ZEROP1 I)
  73.             (COND ((FLOATP R) (SETQ $NUMER T)))
  74.             (SETQ R ($ASKSIGN R))
  75.             (COND ((EQUAL R '$POS) (RETURN 0))
  76.               ((EQUAL R '$NEG) (RETURN (SIMPLIFY '$%PI)))
  77.               (T (MERROR "ATAN(0//0) has been generated."))))
  78.            ((AND (AMONG '%COS R) (AMONG '%SIN I))
  79.             (SETQ VAR 'XZ)
  80.             (NUMDEN (DIV* R I))
  81.             (COND ((AND (EQ (CAAR NN*) '%COS) (EQ (CAAR DN*) '%SIN))
  82.                (RETURN (CADR NN*))))))
  83.          (SETQ A ($SIGN R) B ($SIGN I))
  84.          (COND ((EQ A '$POS) (SETQ A 1))
  85.            ((EQ A '$NEG) (SETQ A -1))
  86.            ((EQ A '$ZERO) (SETQ A 0)))
  87.          (COND ((EQ B '$POS) (SETQ B 1))
  88.            ((EQ B '$NEG) (SETQ B -1))
  89.            ((EQ A '$ZERO) (SETQ B 0)))
  90.          (COND ((EQUAL I 0)
  91.             (RETURN (IF (EQUAL A 1) 0 (SIMPLIFY '$%PI))))
  92.            ((EQUAL R 0)
  93.             (RETURN (COND ((EQUAL B 1) (SIMPLIFY HALF%PI))
  94.                   (T (MUL2 '((RAT SIMP) -1 2)
  95.                        (SIMPLIFY '$%PI)))))))
  96.          (SETQ R (SIMPTIMES (LIST '(MTIMES) A B (DIV* I R)) 1 NIL))
  97.          (RETURN (COND ((ONEP1 R)
  98.                 (ARCHK A B (LIST '(MTIMES) '((RAT) 1 4) '$%PI)))
  99.                ((ALIKE1 R '((MEXPT) 3 ((RAT) 1 2)))
  100.                 (ARCHK A B (LIST '(MTIMES) '((RAT) 1 3) '$%PI)))
  101.                ((ALIKE1 R '((MEXPT) 3 ((RAT) -1 2)))
  102.                 (ARCHK A B (LIST '(MTIMES) '((RAT) 1 6) '$%PI))))))))
  103.  
  104. (declare-top (SPLITFILE BINOML))
  105.  
  106. (DEFMFUN SIMPBINOCOEF (X VESTIGIAL Z) 
  107.  VESTIGIAL ;Ignored.
  108.  (TWOARGCHECK X)
  109.  (LET ((U (SIMPCHECK (CADR X) Z))
  110.        (V (SIMPCHECK (CADDR X) Z))
  111.        (Y))
  112.    (COND ((INTEGERP V)
  113.       (COND ((MINUSP V)
  114.          (IF (AND (INTEGERP U) (MINUSP U) (LESSP V U)) (BINCOMP U (*DIF U V)) 0))
  115.         ((OR (ZEROP V) (EQUAL U V)) 1)
  116.         ((AND (INTEGERP U) (NOT (MINUSP U))) (BINCOMP U (MIN V (*DIF U V))))
  117.         (T (BINCOMP U V))))
  118.      ((INTEGERP (SETQ Y (SUB U V))) (BINCOMP U Y))
  119.      ((AND (FLOATP U) (FLOATP V)) ($MAKEGAMMA (LIST '(%BINOMIAL) U V)))
  120.      (T (EQTEST (LIST '(%BINOMIAL) U V) X)))))
  121.  
  122. (DEFUN BINCOMP (U V) 
  123.        (COND ((MINUSP V) 0)
  124.          ((ZEROP V) 1)
  125.          ((MNUMP U) (BINOCOMP U V))
  126.          (T (MULN (BINCOMP1 U V) NIL)))) 
  127.  
  128. (DEFUN BINCOMP1 (U V) 
  129.  (IF (EQUAL V 1)
  130.      (NCONS U)
  131.      (LIST* U (LIST '(MEXPT) V -1) (BINCOMP1 (ADD2 -1 U) (SUB1 V)))))
  132.  
  133. (DEFMFUN BINOCOMP (U V) 
  134.        (PROG (ANS) 
  135.          (SETQ ANS 1)
  136.     LOOP (IF (ZEROP V) (RETURN ANS))
  137.          (SETQ ANS (TIMESK (TIMESK U ANS) (SIMPLIFY (LIST '(RAT) 1 V))))
  138.          (SETQ U (ADDK -1 U) V (SUB1 V))
  139.          (GO LOOP)))
  140.  
  141. (declare-top (SPLITFILE GAMMA) (SPECIAL $NUMER $GAMMALIM))
  142.  
  143. (DEFMVAR $BETA_ARGS_SUM_TO_INTEGER NIL)
  144.  
  145. (DEFMFUN SIMPBETA (X VESTIGIAL Z &AUX CHECK)
  146.  VESTIGIAL ;Ignored.
  147.  (TWOARGCHECK X)
  148.  (SETQ CHECK X)
  149.  (LET ((U (SIMPCHECK (CADR X) Z)) (V (SIMPCHECK (CADDR X) Z)))
  150.       (COND ((OR (ZEROP1 U) (ZEROP1 V))
  151.          (IF ERRORSW (THROW 'ERRORSW T) (MERROR "Zero argument to BETA")))
  152.         ((OR (AND (FLOATP U) (FLOATP V))
  153.          (AND $NUMER (NUMBERP U) (NUMBERP V)))
  154.          ($MAKEGAMMA (LIST '($BETA) U V)))
  155.         ((OR (AND (INTEGERP U) (PLUSP U)) (AND (INTEGERP V) (PLUSP V)))
  156.          (SETQ X (ADD2 U V))
  157.          (POWER (MUL2 (SUB X 1)
  158.               (SIMPLIFYA (LIST '(%BINOMIAL)
  159.                        (SUB X 2)
  160.                        (SUB (IF (AND (INTEGERP U) (PLUSP U)) U V) 1))
  161.                      T))
  162.             -1))
  163.         ((AND (INTEGERP U) (INTEGERP V))
  164.          (MUL2* (DIV* (LIST '(MFACTORIAL) (SUB1 U))
  165.               (LIST '(MFACTORIAL) (PLUS U V -1)))
  166.             (LIST '(MFACTORIAL) (SUB1 V))))
  167.         ((OR (AND (RATNUMP U) (RATNUMP V) (INTEGERP (SETQ X (ADDK U V))))
  168.          (AND $BETA_ARGS_SUM_TO_INTEGER
  169.               (INTEGERP (SETQ X (EXPAND1 (ADD2 U V) 1 1)))))
  170.          (LET ((W (IF (SYMBOLP V) V U)))
  171.           (DIV* (MUL2* '$%PI
  172.                    (LIST '(%BINOMIAL)
  173.                      (ADD2 (SUB1 X) (NEG W))
  174.                      (SUB1 X)))
  175.             `((%SIN) ((MTIMES) ,W $%PI)))))
  176.         (T (EQTEST (LIST '($BETA) U V) CHECK)))))
  177.  
  178. (DEFMFUN SIMPGAMMA (X VESTIGIAL Z)
  179.  VESTIGIAL ;Ignored.
  180.  (ONEARGCHECK X)
  181.  (LET ((J (SIMPCHECK (CADR X) Z)))
  182.    (COND ((FLOATP J) (GAMMAFLOAT J))
  183.      ((OR (NOT (MNUMP J))
  184.           (RATGREATERP (SIMPABS (LIST '(%ABS) J) 1 T) $GAMMALIM))
  185.       (EQTEST (LIST '(%GAMMA) J) X))
  186.      ((INTEGERP J)
  187.       (COND ((GREATERP J 0) (SIMPFACT (LIST '(MFACTORIAL) (SUB1 J)) 1 NIL))
  188.         (ERRORSW (THROW 'ERRORSW T))
  189.         (T (MERROR "GAMMA(~:M) is undefined" J))))
  190.      ($NUMER (GAMMAFLOAT (FPCOFRAT J)))
  191.      ((ALIKE1 J '((RAT) 1 2))
  192.       (LIST '(MEXPT SIMP) '$%PI J))
  193.      ((OR (RATGREATERP J 1) (RATGREATERP 0 J)) (GAMMARED J))
  194.      (T (EQTEST (LIST '(%GAMMA) J) X)))))
  195.  
  196. (declare-top (flonum sum))
  197.  
  198. (defun gamma (y) ;;; numerical evaluation for 0 < y < 1
  199.        (prog (sum coefs)
  200.          (setq coefs '(0.035868343 -0.193527817
  201.                        0.48219939
  202.                        -0.75670407
  203.                        0.91820685
  204.                        -0.89705693
  205.                        0.98820588
  206.                        -0.57719165))
  207.          (or (atom y) (setq y (fpcofrat y)))
  208.          (setq sum (car coefs) coefs (cdr coefs))
  209.     loop (setq sum (+$ (*$ sum y) (car coefs)))
  210.          (and (setq coefs (cdr coefs)) (go loop))
  211.          (return (+$ (//$ 1.0 y) sum))))
  212.  
  213. (declare-top (notype sum))
  214.  
  215. (defun gammared (a)                    ;A is assumed to
  216.        (prog (m q n)                    ;be '((RAT) M N)
  217.          (cond ((floatp a) (return (gammafloat a))))
  218.          (setq m (cadr a)                ;Numerator
  219.            n (caddr a)                ;denominator
  220.            q (abs (*quo m n)))            ;integer part
  221.          (cond ((minusp m)
  222.             (setq q (add1 q) m (plus m (times n q)))
  223.             (return
  224.              (simptimes (list '(mtimes)
  225.                       (list '(mexpt) n q)
  226.                       (simpgamma (list '(%gamma)
  227.                                (list '(rat) m n))
  228.                          1.
  229.                          nil)
  230.                       (list '(mexpt) (gammac m n q) -1.))
  231.                 1.
  232.                 nil))))
  233.          (return (m* (gammac m n q)
  234.              (simpgamma (list '(%gamma)
  235.                       (list '(rat) (remainder m n) n))
  236.                     1 nil)
  237.              (m^ n (minus q))))))
  238.  
  239. (defun gammac (m n q)
  240.        (do ((ans 1))
  241.        ((lessp q 1) ans)
  242.        (setq q (sub1 q) m (*dif m n) ans (times m ans))))
  243.  
  244. #+nil
  245. (declare-top (flonum a r))
  246.  
  247. #+nil
  248. (defun gammafloat (a) 
  249.        (cond ((= a 1.0) 1.0)
  250.          ((= a 0.0) (merror "GAMMA(0.0) has been generated."))    
  251.          ((and (> a 0.0) (> 1.0 a)) (gamma a))
  252.          ((or (> a 34.82) (< a -34.12))
  253.           (merror "GAMMA(~A) - arithmetic overflow" a))
  254.          (t (do-gammafloat a))))
  255.  
  256. ;
  257. ;(defun gammafloat (a)
  258. ;  (cond ((= a 1.0) 1.0)
  259. ;    ((= a 0.0) (merror "GAMMA(0.0) has been generated."))    
  260. ;    ((and (> a 0.0) (> 1.0 a)) (gamma a))
  261. ;    (t (condition-case ()
  262. ;           (do-gammafloat a)
  263. ;         (dbg:floating-exponent-overflow 
  264. ;           (merror "GAMMA(~A) - arithmetic overflow" a))))))
  265.  
  266. #+nil
  267. (defun do-gammafloat (a)
  268.   (do ((r 1.0 (*$ z r))
  269.        (s (minusp a)) (z (abs a)))
  270.       ((not (greaterp z 1.0))
  271.        (setq r (*$ r (gamma z)))
  272.        (cond (s (t//$ -3.141592654 (*$ a r (sin (*$ 3.141592654 a))) 'gamma))
  273.          (t r))) 
  274.     (setq z (1-$ z))))
  275.  
  276. ;; This implementation is based on Lanczos convergent formula for the
  277. ;; gamma function for Re(z) > 0.  We can use the reflection formula
  278. ;;
  279. ;;    -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
  280. ;;
  281. ;; to handle the case of Re(z) <= 0.
  282. ;;
  283. ;; See http://winnie.fit.edu/~gabdo/gamma.m for some matlab code to
  284. ;; compute this and http://winnie.fit.edu/~gabdo/gamma.txt for a nice
  285. ;; discussion of Lanczos method and an improvement of Lanczos method.
  286. ;;
  287. ;;
  288. ;; The document says this should give about 15 digits of accuracy for
  289. ;; double-precision IEEE floats.  The document also indicates how to
  290. ;; compute a new set of coefficients if you need more range or
  291. ;; accuracy.
  292.  
  293. (defun gamma-lanczos (z)
  294.   (declare (type (complex double-float) z)
  295.        (optimize (safety 3)))
  296.   (let ((g 607/128)
  297.     (c (make-array 15 :element-type 'double-float
  298.                :initial-contents
  299.                '(0.99999999999999709182d0
  300.              57.156235665862923517d0
  301.              -59.597960355475491248d0
  302.              14.136097974741747174d0
  303.              -0.49191381609762019978d0
  304.              .33994649984811888699d-4
  305.              .46523628927048575665d-4
  306.              -.98374475304879564677d-4
  307.              .15808870322491248884d-3
  308.              -.21026444172410488319d-3
  309.              .21743961811521264320d-3
  310.              -.16431810653676389022d-3
  311.              .84418223983852743293d-4
  312.              -.26190838401581408670d-4
  313.              .36899182659531622704d-5))))
  314.     (declare (type (rational 607/128 607/128) g)
  315.          (type (simple-array double-float (15)) c))
  316.     (if (minusp (realpart z))
  317.     ;; Use the reflection formula
  318.     ;; -z*Gamma(z)*Gamma(-z) = pi/sin(pi*z)
  319.     ;; or
  320.     ;; Gamma(z) = pi/z/sin(pi*z)/Gamma(-z)
  321.     ;;
  322.     ;; If z is a negative integer, Gamma(z) is infinity.  Should
  323.     ;; we test for this?  Throw an error?  What
  324.     (/ (float pi 1d0)
  325.        (* (- z) (sin (* (float pi 1d0)
  326.                 z))
  327.           (gamma-lanczos (- z))))
  328.     (let* ((z (- z 1))
  329.            (zh (+ z 1/2))
  330.            (zgh (+ zh 607/128))
  331.            (zp (expt zgh (/ zh 2)))
  332.            (ss 
  333.          (do ((sum 0d0)
  334.               (pp (1- (length c)) (1- pp)))
  335.              ((< pp 1)
  336.               sum)
  337.            (incf sum (/ (aref c pp) (+ z pp))))
  338.          ))
  339.       (* (sqrt (float (* 2 pi) 1d0))
  340.          (+ ss (aref c 0))
  341.          (* zp (exp (- zgh)) zp))))))
  342.  
  343. (defun gammafloat (a)
  344.   (realpart (gamma-lanczos (complex a 0d0))))
  345.  
  346. (declare-top (notype a r))
  347.  
  348. (declare-top (SPLITFILE ERF) (SPECIAL $NUMER $TRIGSIGN))
  349.  
  350. (defmfun simperf (x vestigial z &aux y)
  351.  vestigial ;Ignored.
  352.  (oneargcheck x)
  353.  (setq y (simpcheck (cadr x) z))
  354.  (cond ((zerop1 y) y)
  355.        ((or (floatp y) (and $numer (integerp y))) (erf (float y)))
  356.        ((eq y '$inf) 1)
  357.        ((eq y '$minf) -1)
  358.        ((and $trigsign (mminusp* y)) (neg (list '(%erf simp) (neg y))))
  359.        (t (eqtest (list '(%erf) y) x))))
  360.  
  361. #+nil
  362. (defmfun erf (y)
  363.        (cond ((> (abs y) 4.0) (cond ((> y 0.0) 1.0) (t -1.0)))
  364.          (t ((lambda (t1 xf)
  365.            (declare (flonum t1 xf))
  366.            (setq t1 (//$ (1+$ (*$ xf 0.3275911))))
  367.            (setq 
  368.             t1
  369.             (-$
  370.              1.0
  371.              (*$
  372.               (exp (minus (*$ xf xf)))
  373.               (*$ (+$ (*$ (+$ (*$ (+$ (*$ (+$ (*$ t1
  374.                               1.06140543)
  375.                               -1.45315203)
  376.                           t1)
  377.                           1.42141373)
  378.                       t1)
  379.                       -0.28449674)
  380.                   t1)
  381.                   0.25482959)
  382.               t1))))
  383.            (cond ((> y 0.0) t1) (t (minus t1))))
  384.           0.0
  385.           (abs y)))))
  386.  
  387. (defmfun erf (y)
  388.   (slatec:derf (float y 1d0)))
  389. (defmfun erfc (y)
  390.   (slatec:derfc (float y 1d0)))
  391.  
  392.  
  393. (declare-top (SPLITFILE EMATRIX))
  394.  
  395. (DEFMFUN $ZEROMATRIX (M N) ($EMATRIX M N 0 1 1))
  396.  
  397. (DEFMFUN $EMATRIX (M N VAR I J)
  398.        (PROG (ANS ROW) 
  399.        (COND ((EQUAL M 0) (RETURN (NCONS '($MATRIX SIMP))))
  400.          ((AND (EQUAL N 0) (FIXNUMP M) (> M 0))
  401.           (RETURN (CONS '($MATRIX SIMP) (LIST-OF-MLISTS M))))
  402.          ((NOT (AND (FIXNUMP M) (FIXNUMP N)
  403.             (FIXNUMP I) (FIXNUMP J)
  404.             (> M 0) (> N 0) (> I 0) (> J 0)))
  405.           (MERROR "Incorrect argument to EMATRIX:~%~M"
  406.               (LIST '(MLIST SIMP) M N I J) )))
  407.     LOOP (COND ((= M I) (SETQ ROW (ONEN J N VAR 0)) (GO ON))
  408.            ((ZEROP M) (RETURN (CONS '($MATRIX) (MXC ANS)))))
  409.          (SETQ ROW NIL)
  410.          (DO ((N N (f1- N))) ((ZEROP N)) (SETQ ROW (CONS 0 ROW)))
  411.     ON   (SETQ ANS (CONS ROW ANS) M (f1- M))
  412.          (GO LOOP)))
  413.  
  414. (DEFUN LIST-OF-MLISTS (N)
  415.  (DO ((N N (f1- N)) (L NIL (CONS (NCONS '(MLIST SIMP)) L))) ((= N 0) L)))
  416.  
  417. (declare-top (SPLITFILE COEFM) (SPECIAL $RATMX))
  418.  
  419. (DEFMFUN $COEFMATRIX (EQL VARL) (COEFMATRIX EQL VARL NIL))
  420.  
  421. (DEFMFUN $AUGCOEFMATRIX (EQL VARL) (COEFMATRIX EQL VARL T))
  422.  
  423. (DEFUN COEFMATRIX (EQL VARL IND)
  424.        (PROG (ANS ROW A B ELEM)
  425.          (IF (NOT ($LISTP EQL)) (IMPROPER-ARG-ERR EQL '$COEFMATRIX))
  426.          (IF (NOT ($LISTP VARL)) (IMPROPER-ARG-ERR VARL '$COEFMATRIX))
  427.          (DOLIST (V (CDR VARL))
  428.            (IF (AND (NOT (ATOM V)) (MEMQ (CAAR V) '(MPLUS MTIMES)))
  429.            (MERROR "Improper variable to COEFMATRIX:~%~M" V)))
  430.          (SETQ EQL (NREVERSE (MAPCAR #'MEQHK (CDR EQL)))
  431.            VARL (REVERSE (CDR VARL)))
  432.     LOOP1(IF (NULL EQL) (RETURN (CONS '($MATRIX) (MXC ANS))))
  433.          (SETQ A (CAR EQL) EQL (CDR EQL) ROW NIL)
  434.          (IF IND (SETQ ROW (CONS (CONST1 A VARL) ROW)))
  435.          (SETQ B VARL)
  436.     LOOP2(SETQ ELEM (RATCOEF A (CAR B)))
  437.          (SETQ ROW (CONS (IF $RATMX ELEM (RATDISREP ELEM)) ROW))
  438.          (IF (SETQ B (CDR B)) (GO LOOP2))
  439.          (SETQ ANS (CONS ROW ANS))
  440.          (GO LOOP1)))
  441.  
  442. (DEFUN CONST1 (E VARL) (DOLIST (V VARL) (SETQ E (MAXIMA-SUBSTITUTE 0 V E))) E)
  443.  
  444.  
  445. (declare-top (SPLITFILE ENTERM))
  446.  
  447. (DEFMFUN $ENTERMATRIX (ROWS COLUMNS)
  448.        (PROG (ROW COLUMN VECTOR MATRIX SYM SYMVECTOR)
  449.          (COND ((OR (NOT (FIXNUMP ROWS))
  450.             (NOT (FIXNUMP COLUMNS)))
  451.             (MERROR "ENTERMATRIX called with non-integer arguments")))
  452.          (SETQ ROW 0)
  453.          (COND ((NOT (= ROWS COLUMNS)) (SETQ SYM NIL) (GO OLOOP)))
  454.     QUEST(PRINC "
  455. Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. General
  456. Answer 1, 2, 3 or 4 : ")         (SETQ SYM (RETRIEVE NIL NIL))
  457.          (COND ((NOT (zl-MEMBER SYM '(1 2 3 4))) (GO QUEST)))
  458.     OLOOP(COND ((> (SETQ ROW (f1+ ROW)) ROWS)
  459.             (format t "~%Matrix entered.~%")
  460.             (RETURN (CONS '($MATRIX) (MXC MATRIX)))))
  461.          (COND ((EQUAL SYM 1)
  462.             (SETQ COLUMN ROW)
  463.             (PRINC "Row ") (PRINC ROW) (PRINC " Column ")
  464.             (PRINC COLUMN) (PRINC ":  ") 
  465.             (SETQ MATRIX
  466.              (NCONC MATRIX
  467.               (NCONS (ONEN ROW COLUMNS (MEVAL (RETRIEVE NIL NIL)) 0))))
  468.             (GO OLOOP))
  469.            ((EQUAL SYM 2)
  470.             (SETQ COLUMN (f1- ROW))
  471.             (COND ((EQUAL ROW 1) (GO ILOOP)))
  472.             (SETQ SYMVECTOR 
  473.                    (CONS (NTHCDR COLUMN VECTOR) SYMVECTOR)
  474.                   VECTOR (NREVERSE (MAPCAR 'CAR SYMVECTOR))
  475.               SYMVECTOR (MAPCAR 'CDR SYMVECTOR))
  476.             (GO ILOOP))
  477.            ((EQUAL SYM 3)
  478.             (SETQ COLUMN ROW)
  479.             (COND ((EQUAL ROW 1) (SETQ VECTOR (NCONS 0)) (GO ILOOP)))
  480.             (SETQ SYMVECTOR
  481.               (CONS (MAPCAR 'NEG
  482.                     (NTHCDR (f1- COLUMN) VECTOR))
  483.                 SYMVECTOR)
  484.               VECTOR (NRECONC (MAPCAR 'CAR SYMVECTOR) (NCONS 0))
  485.               SYMVECTOR (MAPCAR 'CDR SYMVECTOR))
  486.             (GO ILOOP)))         
  487.          (SETQ COLUMN 0 VECTOR NIL)
  488.     ILOOP(COND ((> (SETQ COLUMN (f1+ COLUMN)) COLUMNS)
  489.             (SETQ MATRIX (NCONC MATRIX (NCONS VECTOR)))
  490.             (GO OLOOP)))
  491.          (PRINC "Row ") (PRINC ROW) (PRINC " Column ")
  492.          (PRINC COLUMN) (PRINC ":  ") 
  493.          (SETQ VECTOR (NCONC VECTOR (NCONS (MEVAL (RETRIEVE NIL NIL)))))
  494.          (GO ILOOP)))
  495.  
  496. (declare-top (splitfile xthru) (special sn* sd* rsn*))
  497.  
  498. (DEFMFUN $xthru (e)
  499.        (cond ((atom e) e)
  500.          ((mtimesp e) (muln (mapcar '$xthru (cdr e)) nil))
  501.          ((mplusp e) (simplify (comdenom (mapcar '$xthru (cdr e)) t)))
  502.          ((mexptp e) (power ($xthru (cadr e)) (caddr e)))
  503.          ((memq (caar e) '(mequal mlist $matrix))
  504.           (cons (car e) (mapcar '$xthru (cdr e))))
  505.          (t e))) 
  506.  
  507. (defun comdenom (l ind) 
  508.   (prog (n d) 
  509.     (prodnumden (car l))
  510.     (setq n (m*l sn*) sn* nil)
  511.     (setq d (m*l sd*) sd* nil)
  512.    loop    (setq l (cdr l))
  513.     (cond ((null l)
  514.            (return (cond (ind (div* (cond (rsn* ($ratsimp n))
  515.                           (t n))
  516.                     d))
  517.                  (t (list n d))))))
  518.     (prodnumden (car l))
  519.     (setq d (comdenom1 n d (m*l sn*) (m*l sd*)))
  520.     (setq n (car d))
  521.     (setq d (cadr d))
  522.     (go loop)))
  523.  
  524. (defun prodnumden (e) 
  525.  (cond ((atom e) (prodnd (list e)))
  526.        ((eq (caar e) 'mtimes) (prodnd (cdr e)))
  527.        (t (prodnd (list e)))))
  528.  
  529. (defun prodnd (l) 
  530.      (prog (e) 
  531.            (setq l (reverse l))
  532.            (setq sn* nil sd* nil)
  533.       loop (cond ((null l) (return nil)))
  534.            (setq e (car l))
  535.            (cond ((atom e) (setq sn* (cons e sn*)))
  536.              ((ratnump e)
  537.               (cond ((not (equal 1. (cadr e)))
  538.                  (setq sn* (cons (cadr e) sn*))))
  539.               (setq sd* (cons (caddr e) sd*)))
  540.              ((and (eq (caar e) 'mexpt)
  541.                (mnegp (caddr e)))
  542.               (setq sd* (cons (power (cadr e)
  543.                          (timesk -1 (caddr e)))
  544.                       sd*)))
  545.              (t (setq sn* (cons e sn*))))
  546.            (setq l (cdr l))
  547.            (go loop)))
  548.  
  549. (defun comdenom1 (a b c d) 
  550.        (prog (b1 c1) 
  551.          (prodnumden (div* b d))
  552.          (setq b1 (m*l sn*) sn* nil)
  553.          (setq c1 (m*l sd*) sd* nil)
  554.          (return
  555.           (list (add2 (m* a c1) (m* c b1))
  556.             (mul2 d b1)))))
  557.  
  558. (declare-top (SPLITFILE XRTOUT)
  559.      (SPECIAL $GLOBALSOLVE $BACKSUBST $DISPFLAG $NOLABELS
  560.           $LINSOLVE_PARAMS $%RNUM_LIST AX LINELABLE $LINECHAR 
  561.           $LINENUM SOL *MOSESFLAG) 
  562.      (FIXNUM TIM $LINENUM))
  563.  
  564. (DEFUN XRUTOUT (AX N M VARL IND)
  565.  (LET (($LINSOLVE_PARAMS (AND $BACKSUBST $LINSOLVE_PARAMS)))
  566.   (PROG (IX IMIN J ANS ZZ M-1 SOL TIM CHK ZZZ)
  567.     (SETQ AX (GET-ARRAY-POINTER AX) TIM 0)
  568.     (IF $LINSOLVE_PARAMS (SETQ $%RNUM_LIST (LIST '(MLIST))))
  569.     (SETQ IMIN (MIN (SETQ M-1 (f1- M)) N))
  570.     (SETQ IX (MAX IMIN (LENGTH VARL)))
  571.    LOOP (IF (ZEROP IX) (IF IND (GO OUT) (RETURN (CONS '(MLIST) ZZ))))
  572.         (WHEN (OR (> IX IMIN) (EQUAL (CAR (ARRAYCALL T AX IX IX)) 0))
  573.           (STORE (ARRAYCALL T AX 0 IX)
  574.              (RFORM (IF $LINSOLVE_PARAMS (MAKE-PARAM) (ITH VARL IX))))
  575.           (IF $LINSOLVE_PARAMS (GO SAVAL) (GO NEXT)))
  576.     (SETQ ANS (ARRAYCALL T AX IX M))
  577.     (STORE (ARRAYCALL T AX IX M) NIL)
  578.     (DO ((J (f1+ IX) (f1+ J))) ((> J M-1))
  579.         (SETQ ANS (RATDIF ANS (RATTIMES (ARRAYCALL T AX IX J) 
  580.                         (ARRAYCALL T AX 0 J)
  581.                         T)))
  582.         (STORE (ARRAYCALL T AX IX J ) NIL))
  583.     (STORE (ARRAYCALL T AX 0 IX) (RATQUOTIENT ANS (ARRAYCALL T AX IX IX)))
  584.     (STORE (ARRAYCALL T AX IX IX ) NIL)
  585.     (SETQ ANS NIL)
  586.    SAVAL(PUSH (COND (*MOSESFLAG (ARRAYCALL T AX 0 IX)) 
  587.             (T (LIST (IF $GLOBALSOLVE '(MSETQ) '(MEQUAL))
  588.                  (ITH VARL IX)
  589.                  (SIMPLIFY (RDIS (ARRAYCALL T AX 0 IX))))))
  590.           ZZ)
  591.         (IF (NOT $BACKSUBST)
  592.         (STORE (ARRAYCALL T AX 0 IX) (RFORM (ITH VARL IX))))
  593.     (AND $GLOBALSOLVE (MEVAL (CAR ZZ)))
  594.    NEXT (SETQ IX (f1- IX))
  595.     (GO LOOP)
  596.    OUT
  597.     (COND ($DISPFLAG (MTELL "Solution~%")))
  598.     (SETQ J 1 SOL (LIST '(MLIST)) CHK (CHECKLABEL $LINECHAR))
  599.     (DO ((LL ZZ (CDR LL))) ((NULL LL)) (SETQ ZZZ (CAR LL))
  600.        (SETQ ZZZ (LIST '(MLABLE)
  601.                (PROGN (IF CHK (SETQ CHK NIL)
  602.                       (SETQ $LINENUM (f1+ $LINENUM)))
  603.                   ((LAMBDA ($NOLABELS)(MAKELABEL $LINECHAR))
  604.                    (AND $NOLABELS $GLOBALSOLVE))
  605.                   LINELABLE)
  606.                (COND ((NOT (AND $NOLABELS $GLOBALSOLVE))
  607.                   (SET LINELABLE ZZZ)))))
  608.        (NCONC SOL (NCONS LINELABLE))
  609.        (COND ($DISPFLAG (SETQ TIM (RUNTIME))
  610.                 (MTELL-OPEN "~%~M" ZZZ)
  611.                 (TIMEORG TIM))
  612.          (T (PUTPROP LINELABLE T 'NODISP))))
  613.     (RETURN SOL))))
  614.