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 / asum.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  24.5 KB  |  765 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 asum)
  13. (load-macsyma-macros rzmac)
  14.  
  15. (declare-top(special opers *a *n $factlim sum msump *i
  16.           *opers-list opers-list $ratsimpexpons makef)
  17.      (*expr sum)
  18.      (fixnum %n %k %i %m $genindex)
  19.      (genprefix sm))
  20.  
  21. (sloop for (x y) on 
  22.        '(%cot %tan %csc %sin %sec %cos %coth %tanh %csch %sinh %sech %cosh)
  23.        by 'cddr do (putprop x y 'recip) (putprop y x 'recip))
  24.  
  25. (defun nill () '(nil))
  26.  
  27. (defmvar $zeta%pi t)
  28.  
  29.     (comment Polynomial predicates and other such things)
  30.  
  31. (defun poly? (exp var)
  32.     (cond ((or (atom exp) (free exp var)))
  33.       ((memq (caar exp) '(mtimes mplus))
  34.        (do ((exp (cdr exp) (cdr exp)))
  35.            ((null exp) t)
  36.            (and (null (poly? (car exp) var)) (return nil))))
  37.       ((and (eq (caar exp) 'mexpt)
  38.          (integerp (caddr exp))
  39.         (> (caddr exp) 0))
  40.        (poly? (cadr exp) var))))
  41.  
  42. (defun smono (x var) (smonogen x var t))
  43.  
  44. (defun smonop (x var) (smonogen x var nil))
  45.  
  46. (defun smonogen (x var fl)  ; fl indicates whether to return *a *n
  47.     (cond ((free x var) (and fl (setq *n 0 *a x)) t)
  48.       ((atom x) (and fl (setq *n (setq *a 1))) t)
  49.       ((eq (caar x) 'mtimes)
  50.        (do ((x (cdr x) (cdr x))
  51.         (a '(1)) (n '(0)))
  52.            ((null x)
  53.         (and fl (setq *n (addn n nil) *a (muln a nil))) t)
  54.            (let (*a *n)
  55.          (if (smonogen (car x) var fl)
  56.              (and fl (setq a (cons *a a) n (cons *n n)))
  57.              (return nil)))))
  58.       ((eq (caar x) 'mexpt)
  59.        (cond ((and (free (caddr x) var) (eq (cadr x) var))
  60.           (and fl (setq *n (caddr x) *a 1)) t)))))
  61.  
  62.     (comment Factorial Stuff)
  63.  
  64. (setq $factlim -1 makef nil)
  65.  
  66. (DEFMFUN $genfact (&rest l)
  67.   (cons '(%genfact) (copy-rest-arg l)))
  68.  
  69. (defun gfact (n %m i)
  70.   (cond ((minusp %m) (improper-arg-err %m '$genfact))
  71.     ((= %m 0) 1)
  72.     (t (prog (ans)
  73.          (setq ans n)
  74.            a (if (= %m 1) (return ans))
  75.              (setq n (m- n i) %m (f1- %m) ans (m* ans n))
  76.          (go a)))))
  77.  
  78. ;(defun factorial (%i)
  79. ;    (cond ((< %i 2) 1)
  80. ;          (t (do ((x 1 (times %i x)) (%i %i (f1- %i)))
  81. ;             ((= %i 1) x)))))
  82. ;; people like to do big factorials so do them a bit faster.
  83. ;; by breaking into chunks.
  84.  
  85. (defun factorial (n &aux (ans 1)  )
  86.   (let* ((vec (make-array (if (< n 100) 1 20) :initial-element 1))
  87.      (m (length vec))
  88.      (j 0))
  89.   (sloop for i from 1 to n
  90.      do (setq j (mod i m))
  91.      (setf (aref vec j) (* (aref vec j) i)))
  92.   (sloop for v in-array vec
  93.      do (setq ans (* ans v)))
  94.   ans))
  95.  
  96. (defmfun simpfact (x y z)
  97.     (oneargcheck x)
  98.     (setq y (simpcheck (cadr x) z))
  99.     (cond ((or (floatp y) (and (not makef) (ratnump y) (equal (caddr y) 2)))
  100.            (simplifya (makegamma1 (list '(mfactorial) y)) nil))
  101.           ((or (not (fixnump y)) (not (> y -1)))
  102.            (eqtest (list '(mfactorial) y) x))
  103.           ((or (minusp $factlim) (not (greaterp y $factlim)))
  104.            (factorial y))
  105.           (t (eqtest (list '(mfactorial) y) x))))
  106.  
  107. (defun makegamma1 (e)
  108.   (cond ((atom e) e)
  109.     ((eq (caar e) 'mfactorial)
  110.      (list '(%gamma) (list '(mplus) 1 (makegamma1 (cadr e)))))
  111.     ((eq (caar e) '%elliptic_kc)
  112.      ;; Complete elliptic integral of the first kind
  113.      (cond ((alike1 (cadr e) '((rat simp) 1 2))
  114.         ;; K(1/2) = gamma(1/4)/4/sqrt(pi)
  115.         '((mtimes simp) ((rat simp) 1 4)
  116.           ((mexpt simp) $%pi ((rat simp) -1 2))
  117.           ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2)))
  118.            ((or (alike1 (cadr e)
  119.                 '((mtimes simp) ((rat simp) 1 4)
  120.                   ((mplus simp) 2
  121.                    ((mexpt simp) 3 ((rat simp) 1 2)))))
  122.             (alike1 (cadr e)
  123.                 '((mplus simp) ((rat simp) 1 2)
  124.                   ((mtimes simp) ((rat simp) 1 4)
  125.                    ((mexpt simp) 3 ((rat simp) 1 2)))))
  126.             (alike1 (cadr e)
  127.                 ;; 1/(8-4*sqrt(3))
  128.                 '((mexpt simp)
  129.                   ((mplus simp) 8
  130.                    ((mtimes simp) -4
  131.                 ((mexpt simp) 3 ((rat simp) 1 2))))
  132.                   -1)))
  133.             ;; K((2+sqrt(3)/4))
  134.             '((mtimes simp) ((rat simp) 1 4)
  135.               ((mexpt simp) 3 ((rat simp) 1 4))
  136.               ((mexpt simp) $%pi ((rat simp) -1 2))
  137.               ((%gamma simp) ((rat simp) 1 6))
  138.               ((%gamma simp) ((rat simp) 1 3))))
  139.         ((or (alike1 (cadr e)
  140.                  ;; (2-sqrt(3))/4
  141.                  '((mtimes simp) ((rat simp) 1 4)
  142.                    ((mplus simp) 2
  143.                 ((mtimes simp) -1
  144.                  ((mexpt simp) 3 ((rat simp) 1 2))))))
  145.              (alike1 (cadr e)
  146.                  ;; 1/2-sqrt(3)/4
  147.                  '((mplus simp) ((rat simp) 1 2)
  148.                    ((mtimes simp) ((rat simp) -1 4)
  149.                 ((mexpt simp) 3 ((rat simp) 1 2)))))
  150.              (alike (cadr e)
  151.                 ;; 1/(4*sqrt(3)+8)
  152.                 '((mexpt simp)
  153.                   ((mplus simp) 8
  154.                    ((mtimes simp) 4
  155.                 ((mexpt simp) 3 ((rat simp) 1 2))))
  156.                   -1)))
  157.              ;; K((2-sqrt(3))/4)
  158.              '((mtimes simp) ((rat simp) 1 4)
  159.                ((mexpt simp) 3 ((rat simp) -1 4))
  160.                ((mexpt simp) $%pi ((rat simp) -1 2))
  161.                ((%gamma simp) ((rat simp) 1 6))
  162.                ((%gamma simp) ((rat simp) 1 3))))
  163.          ((or
  164.            ;; (3-2*sqrt(2))/(3+2*sqrt(2))
  165.            (alike1 (cadr e)
  166.                '((mtimes simp)
  167.                  ((mplus simp) 3 
  168.                   ((mtimes simp) -2 
  169.                    ((mexpt simp) 2 ((rat simp) 1 2))))
  170.                  ((mexpt simp)
  171.                   ((mplus simp) 3 
  172.                    ((mtimes simp) 2
  173.                 ((mexpt simp) 2 ((rat simp) 1 2)))) -1)))
  174.            ;; 17 - 12*sqrt(2)
  175.            (alike1 (cadr e)
  176.                '((mplus simp) 17 
  177.                  ((mtimes simp) -12
  178.                   ((mexpt simp) 2 ((rat simp) 1 2)))))
  179.            ;;   (2*SQRT(2) - 3)/(2*SQRT(2) + 3)
  180.            (alike1 (cadr e)
  181.                '((mtimes simp) -1
  182.                  ((mplus simp) -3
  183.                   ((mtimes simp) 2 
  184.                    ((mexpt simp) 2 ((rat simp) 1 2))))
  185.                  ((mexpt simp)
  186.                   ((mplus simp) 3
  187.                    ((mtimes simp) 2
  188.                 ((mexpt simp) 2 ((rat simp) 1 2))))
  189.                   -1))))
  190.           '((mtimes simp) ((rat simp) 1 8)
  191.             ((mexpt simp) 2 ((rat simp) -1 2))
  192.             ((mplus simp) 1 ((mexpt simp) 2 ((rat simp) 1 2)))
  193.             ((mexpt simp) $%pi ((rat simp) -1 2))
  194.             ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2)))
  195.          (t
  196.           ;; Give up
  197.           e)))
  198.     ((eq (caar e) '%elliptic_ec)
  199.      ;; Complete elliptic integral of the second kind
  200.      (cond ((alike1 (cadr e) '((rat simp) 1 2))
  201.         ;; 2*E(1/2) - K(1/2) = 2*%pi^(3/2)*gamma(1/4)^(-2)
  202.         '((mplus simp)
  203.           ((mtimes simp) ((mexpt simp) $%pi ((rat simp) 3 2))
  204.            ((mexpt simp) 
  205.             ((%gamma simp irreducible) ((rat simp) 1 4)) -2))
  206.           ((mtimes simp) ((rat simp) 1 8)
  207.            ((mexpt simp) $%pi ((rat simp) -1 2))
  208.            ((mexpt simp) ((%gamma simp) ((rat simp) 1 4)) 2))))
  209.            ((or (alike1 (cadr e)
  210.                 '((mtimes simp) ((rat simp) 1 4)
  211.                   ((mplus simp) 2
  212.                    ((mtimes simp) -1
  213.                 ((mexpt simp) 3 ((rat simp) 1 2))))))
  214.             (alike1 (cadr e)
  215.                 '((mplus simp) ((rat simp) 1 2)
  216.                   ((mtimes simp) ((rat simp) -1 4)
  217.                    ((mexpt simp) 3 ((rat simp) 1 2))))))
  218.         ;; E((2-sqrt(3))/4)
  219.         ;;
  220.         ;; %pi/4/sqrt(3) = K*(E-(sqrt(3)+1)/2/sqrt(3)*K)
  221.         '((mplus simp)
  222.           ((mtimes simp) ((mexpt simp) 3 ((rat simp) -1 4))
  223.            ((mexpt simp) $%pi ((rat simp) 3 2))
  224.            ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1)
  225.            ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1))
  226.           ((mtimes simp) ((rat simp) 1 8)
  227.            ((mexpt simp) 3 ((rat simp) -3 4))
  228.            ((mexpt simp) $%pi ((rat simp) -1 2))
  229.            ((%gamma simp) ((rat simp) 1 6))
  230.            ((%gamma simp) ((rat simp) 1 3)))
  231.           ((mtimes simp) ((rat simp) 1 8)
  232.            ((mexpt simp) 3 ((rat simp) -1 4))
  233.            ((mexpt simp) $%pi ((rat simp) -1 2))
  234.            ((%gamma simp) ((rat simp) 1 6))
  235.            ((%gamma simp) ((rat simp) 1 3)))))
  236.            ((or (alike1 (cadr e)
  237.                 '((mtimes simp) ((rat simp) 1 4)
  238.                   ((mplus simp) 2
  239.                    ((mexpt simp) 3 ((rat simp) 1 2)))))
  240.             (alike1 (cadr e)
  241.                 '((mplus simp) ((rat simp) 1 2)
  242.                   ((mtimes simp) ((rat simp) 1 4)
  243.                    ((mexpt simp) 3 ((rat simp) 1 2))))))
  244.         ;; E((2+sqrt(3))/4)
  245.         ;;
  246.         ;; %pi*sqrt(3)/4 = K1*(E1-(sqrt(3)-1)/2/sqrt(3)*K1)
  247.         '((mplus simp)
  248.           ((mtimes simp) 3 ((mexpt simp) 3 ((rat simp) -3 4))
  249.            ((mexpt simp) $%pi ((rat simp) 3 2))
  250.            ((mexpt simp) ((%gamma simp) ((rat simp) 1 6)) -1)
  251.            ((mexpt simp) ((%gamma simp) ((rat simp) 1 3)) -1))
  252.           ((mtimes simp) ((rat simp) 3 8)
  253.            ((mexpt simp) 3 ((rat simp) -3 4))
  254.            ((mexpt simp) $%pi ((rat simp) -1 2))
  255.            ((%gamma simp) ((rat simp) 1 6))
  256.            ((%gamma simp) ((rat simp) 1 3)))
  257.           ((mtimes simp) ((rat simp) -1 8)
  258.            ((mexpt simp) 3 ((rat simp) -1 4))
  259.            ((mexpt simp) $%pi ((rat simp) -1 2))
  260.            ((%gamma simp) ((rat simp) 1 6))
  261.            ((%gamma simp) ((rat simp) 1 3)))))
  262.            (t
  263.         e)))
  264.     (t (recur-apply #'makegamma1 e))))
  265.  
  266. (defmfun simpgfact (x vestigial z)
  267.   vestigial ;Ignored.
  268.   (if (not (= (length x) 4)) (wna-err '$genfact))
  269.   (setq z (mapcar #'(lambda (q) (simpcheck q z)) (cdr x)))
  270.   (let ((a (car z)) (b ($entier (cadr z))) (c (caddr z)))
  271.     (cond ((and (fixnump b) (> b -1)) (gfact a b c))
  272.       ((integerp b) (merror "Bad second argument to GENFACT: ~:M" b))
  273.       (t (eqtest (list '(%genfact) a
  274.                (if (and (not (atom b))
  275.                     (eq (caar b) '$entier))
  276.                    (cadr b)
  277.                    b)
  278.                c)
  279.               x)))))
  280.  
  281.     (comment SUM begins)
  282.  
  283. (defmvar $cauchysum nil
  284.      "When multiplying together sums with INF as their upper limit, 
  285. causes the Cauchy product to be used rather than the usual product.
  286. In the Cauchy product the index of the inner summation is a function of 
  287. the index of the outer one rather than varying independently."
  288.      modified-commands '$sum)
  289.  
  290. (defmvar $gensumnum 0
  291.      "The numeric suffix used to generate the next variable of
  292. summation.  If it is set to FALSE then the index will consist only of
  293. GENINDEX with no numeric suffix."
  294.      modified-commands '$sum
  295.      setting-predicate #'(lambda (x) (or (null x) (integerp x))))
  296.  
  297. (defmvar $genindex '$i
  298.      "The alphabetic prefix used to generate the next variable of
  299. summation when necessary."
  300.      modified-commands '$sum
  301.      setting-predicate #'symbolp)
  302.  
  303. (defmvar $zerobern t)
  304. (defmvar $simpsum nil)
  305. (defmvar $sumhack nil)
  306. (defmvar $prodhack nil)
  307.  
  308. (defvar *infsumsimp t)
  309.  
  310. ;; These variables should be initialized where they belong.
  311.  
  312. (setq $wtlevel nil $cflength 1
  313.       $weightlevels '((mlist)) *trunclist nil $taylordepth 3
  314.       $maxtaydiff 4 $verbose nil $psexpand nil ps-bmt-disrep t
  315.       silent-taylor-flag nil)
  316.  
  317. (defmacro sum-arg (sum) `(cadr ,sum))
  318.  
  319. (defmacro sum-index (sum) `(caddr ,sum))
  320.  
  321. (defmacro sum-lower (sum) `(cadddr ,sum))
  322.  
  323. (defmacro sum-upper (sum) `(cadr (cdddr ,sum)))
  324.  
  325. (defmspec $sum (l) (setq l (cdr l))
  326.   (if (= (length l) 4)
  327.       (dosum (car l) (cadr l) (meval (caddr l)) (meval (cadddr l)) t)
  328.       (wna-err '$sum)))
  329.  
  330. (defmspec $lsum (l) (setq l (cdr l))
  331.   (or (= (length l) 3) (wna-err '$lsum))
  332.   (let ((form (car l))
  333.     (ind (cadr l))
  334.     (lis (meval (caddr l)))
  335.     (ans 0))
  336.     (or (symbolp ind) (merror "Second argument not a variable ~M" ind))
  337.     (cond (($listp lis)
  338.        (sloop for v in (cdr lis)
  339.           with lind = (cons ind nil)
  340.           for w = (cons v nil)
  341.           do
  342.           (setq ans (add* ans  (mbinding (lind w) (meval form)))))
  343.        ans)
  344.       (t `((%lsum) ,form ,ind ,lis)))))
  345.     
  346. (defmfun simpsum (x y z)
  347.   (let (($ratsimpexpons t)) (setq y (simplifya (sum-arg x) z)))
  348.   (simpsum1 y (sum-index x) (simplifya (sum-lower x) z)
  349.                 (simplifya (sum-upper x) z)))
  350.  
  351. (defun simpsum1 (exp i lo hi)
  352.   (cond ((not (symbolp i)) (merror "Improper index to SUM:~%~M" i))
  353.     ((equal lo hi) (mbinding ((list i) (list hi)) (meval exp)))
  354.     ((and (atom exp)
  355.           (not (eq exp i))
  356.           (getl '%sum '($outative $linear)))
  357.      (freesum exp lo hi 1))
  358.     ((null $simpsum) (list (get '%sum 'msimpind) exp i lo hi))
  359.     (t (simpsum2 exp i lo hi))))
  360.  
  361. (DEFMFUN DOSUM (EXP IND LOW HI SUMP)
  362.  (IF (NOT (SYMBOLP IND))
  363.      (MERROR "Improper index to ~:M:~%~M" (IF SUMP '$SUM '$PRODUCT) IND))
  364.  (UNWIND-PROTECT 
  365.   (PROG (U *I LIND L*I *HL)
  366.     (WHEN (IF SUMP (NOT $SUMHACK) (NOT $PRODHACK))
  367.           (ASSUME (LIST '(MGEQP) IND LOW))
  368.           (IF (NOT (EQ HI '$INF))
  369.           (ASSUME (LIST '(MGEQP) HI IND))))
  370.     (SETQ LIND (cons IND nil))
  371.     (COND ((NOT (FIXNUMP (SETQ *HL (M- HI LOW))))
  372.            (SETQ EXP (MEVALSUMARG EXP IND))
  373.            (RETURN (CONS (IF SUMP '(%SUM) '(%PRODUCT))
  374.                  (LIST EXP IND LOW HI))))
  375.           ((EQUAL *HL -1) (RETURN (IF SUMP 0 1)))
  376.           ((SIGNP L *HL)
  377.            (COND ((AND SUMP $SUMHACK)
  378.               (RETURN
  379.                (M- (DOSUM EXP IND (M+ 1 HI) (M- LOW 1) T))))
  380.              ((AND (NOT SUMP) $PRODHACK)
  381.               (RETURN
  382.                (M// (DOSUM EXP IND (M+ 1 HI) (M- LOW 1) NIL)))))
  383.            (MERROR "Lower bound to ~:M: ~M~%is greater than the upper bound: ~M"
  384.                (IF SUMP '$SUM '$PRODUCT) LOW HI)))
  385.     (SETQ *I LOW L*I (LIST *I) U (IF SUMP 0 1))
  386.     LO  (SETQ U (SIMPLIFYA
  387.          (LIST (IF SUMP '(MPLUS) '(MTIMES))
  388.                (MBINDING (LIND L*I) (MEVAL EXP)) 
  389.                U)
  390.          T))
  391.     (IF (EQUAL *I HI) (RETURN U))
  392.     (SETQ *I (CAR (RPLACA L*I (M+ *I 1))))
  393.     (GO LO))
  394.   (COND ((IF SUMP (NOT $SUMHACK) (NOT $PRODHACK))
  395.      (FORGET (LIST '(MGEQP) IND LOW))
  396.      (IF (NOT (EQ HI '$INF)) (FORGET (LIST '(MGEQP) HI IND)))))))
  397.  
  398. (defmfun do%sum (l op)
  399.     (if (not (= (length l) 4)) (wna-err op))
  400.     (let ((ind (cadr l)))
  401.       (if (mquotep ind) (setq ind (cadr ind)))
  402.       (if (not (symbolp ind))
  403.       (merror "Illegal index to ~:M:~%~M" op ind))
  404.       (list (mevalsumarg (car l) ind)
  405.         ind (meval (caddr l)) (meval (cadddr l)))))
  406.  
  407. (defun mevalsumarg (exp ind)
  408.   (let ((msump t) (lind (list ind)))
  409.     (mbinding (lind lind)
  410.           (resimplify (mevalatoms (if (and (not (atom exp))
  411.                            (get (caar exp)
  412.                             'mevalsumarg-macro))
  413.                       (funcall (get (caar exp)
  414.                             'mevalsumarg-macro)
  415.                            exp)
  416.                       exp))))))
  417.  
  418.     (comment Multiplication of Sums)
  419.  
  420. (defun gensumindex nil
  421.   (implode (nconc (exploden $genindex)
  422.           (and $gensumnum
  423.                (mexploden (setq $gensumnum (f1+ $gensumnum)))))))
  424.  
  425. (defun sumtimes (x y)
  426.   (cond ((null x) y)
  427.     ((null y) x)
  428.     ((or (safe-zerop  x) (safe-zerop y)) 0)
  429.     ((or (atom x) (not (eq (caar x) '%sum))) (sumultin x y))
  430.     ((or (atom y) (not (eq (caar y) '%sum))) (sumultin y x))
  431.     (t (let (u v i j)
  432.          (if (great (sum-arg x) (sum-arg y)) (setq u y v x) (setq u x v y))
  433.          (setq i (let ((ind (gensumindex)))
  434.                (setq u (subst ind (sum-index u) u)) ind))
  435.          (setq j (let ((ind (gensumindex)))
  436.                (setq v (subst ind (sum-index v) v)) ind))
  437.          (if (and $cauchysum (eq (sum-upper u) '$inf)
  438.                  (eq (sum-upper v) '$inf))
  439.          (list '(%sum)
  440.                (list '(%sum)
  441.                  (sumtimes (MAXIMA-SUBSTITUTE j i (sum-arg u))
  442.                        (MAXIMA-SUBSTITUTE (m- i j) j (sum-arg v)))
  443.                  j (sum-lower u) (m- i (sum-lower v)))
  444.                i (m+ (sum-lower u) (sum-lower v)) '$inf)
  445.          (list '(%sum)
  446.                (list '(%sum) (sumtimes (sum-arg u) (sum-arg v))
  447.                      j (sum-lower v) (sum-upper v))
  448.                i (sum-lower u) (sum-upper u)))))))
  449.  
  450. (defun sumultin (x s)  ; Multiplies x into a sum adjusting indices.
  451.   (cond ((or (atom s) (not (eq (caar s) '%sum))) (m* x s))
  452.     ((free x (sum-index s))
  453.      (list* (car s) (sumultin x (sum-arg s)) (cddr s)))
  454.     (t (let ((ind (gensumindex)))
  455.          (list* (car s)
  456.             (sumultin x (subst ind (sum-index s) (sum-arg s)))
  457.             ind
  458.             (cdddr s))))))
  459.  
  460.     (comment Addition of Sums)
  461.  
  462. (defun sumpls (sum out)
  463.   (prog (l)
  464.     (if (null out) (return (cons sum nil)))
  465.     (setq out (setq l (cons nil out)))
  466.    a     (if (null (cdr out)) (return (cons sum (cdr l))))
  467.     (and (not (atom (cadr out)))
  468.          (consp (caadr out))
  469.          (eq (caar (cadr out)) '%sum)
  470.          (alike1 (sum-arg (cadr out)) (sum-arg sum))
  471.          (alike1 (sum-index (cadr out)) (sum-index sum))
  472.          (cond ((onediff (sum-upper (cadr out)) (sum-lower sum))
  473.             (setq sum (list (car sum)
  474.                     (sum-arg sum)
  475.                     (sum-index sum)
  476.                     (sum-lower (cadr out))
  477.                     (sum-upper sum)))
  478.             (rplacd out (cddr out))
  479.             (go a))
  480.            ((onediff (sum-upper sum) (sum-lower (cadr out)))
  481.             (setq sum (list (car sum)
  482.                     (sum-arg sum)
  483.                     (sum-index sum)
  484.                     (sum-lower sum)
  485.                     (sum-upper (cadr out))))
  486.             (rplacd out (cddr out))
  487.             (go a))))
  488.     (setq out (cdr out))
  489.     (go a)))
  490.  
  491. (defun onediff (x y) (equal 1 (m- y x)))
  492.  
  493. (defun freesum (e b a q) (m* e q (m- (m+ a 1) b)))
  494.  
  495.     (comment Linear operator stuff)
  496.  
  497. (defparameter *opers-list '(($linear . linearize1)))
  498. (defparameter  opers (list '$linear))
  499.  
  500. (defun oper-apply (e z)
  501.        (cond ((null opers-list)
  502.           (let ((w (get (caar e) 'operators)))
  503.         (if w (funcall w e 1 z) (simpargs e z))))
  504.          ((get (caar e) (caar opers-list))
  505.           (let ((opers-list (cdr opers-list))
  506.             (fun (cdar opers-list)))
  507.         (funcall fun e z)))
  508.          (t (let ((opers-list (cdr opers-list)))
  509.           (oper-apply e z)))))
  510.  
  511. (defun linearize1 (e z)  ; z = t means args already simplified.
  512.        (linearize2
  513.     (cons (car e) (mapcar #'(lambda (q) (simpcheck q z)) (cdr e)))
  514.     nil))
  515.  
  516. (defun opident (op) (cond ((eq op 'mplus) 0) ((eq op 'mtimes) 1)))
  517.  
  518. (defun rem-const (e)    ;removes constantp stuff
  519.        (do ((l (cdr e) (cdr l))
  520.         (a (list (opident (caar e))))
  521.         (b (list (opident (caar e)))))
  522.        ((null l)
  523.         (cons (simplifya (cons (list (caar e)) a) nil)
  524.           (simplifya (cons (list (caar e)) b) nil)))
  525.        (if (or (mnump (car l)) (constant (car l)))
  526.            (setq a (cons (car l) a))
  527.            (setq b (cons (car l) b)))))
  528.  
  529. (defun linearize2 (e times)
  530.        (cond ((linearconst e))
  531.          ((atom (cadr e)) (oper-apply e t))
  532.          ((eq (caar (cadr e)) 'mplus)
  533.           (addn (mapcar #'(lambda (q)
  534.                       (linearize2 (list* (car e) q (cddr e)) nil))
  535.                 (cdr (cadr e)))
  536.             t))
  537.          ((and (eq (caar (cadr e)) 'mtimes) (null times))
  538.           (let ((z (if (and (cddr e)
  539.                 (or (atom (caddr e))
  540.                     ($subvarp (caddr e))))
  541.                 (partition (cadr e) (caddr e) 1)
  542.                (rem-const (cadr e))))
  543.             (w))
  544.         (setq w (linearize2 (list* (car e)
  545.                        (simplifya (cdr z) t)
  546.                        (cddr e))
  547.                     t))
  548.         (linearize3 w e (car z))))
  549.          (t (oper-apply e t))))
  550.  
  551. (defun linearconst (e)
  552.   (if (or (mnump (cadr e))
  553.       (constant (cadr e))
  554.       (and (cddr e)
  555.            (or (atom (caddr e)) (memq 'array (cdar (caddr e))))
  556.            (free (cadr e) (caddr e))))
  557.       (if (or (zerop1 (cadr e))
  558.           (and (memq (caar e) '(%sum %integrate))
  559.            (= (length e) 5)
  560.            (or (eq (cadddr e) '$minf)
  561.                (memq (car (cddddr e)) '($inf $infinity)))
  562.            (eq ($asksign (cadr e)) '$zero)))
  563.       0
  564.       (let ((w (oper-apply (list* (car e) 1 (cddr e)) t)))
  565.         (linearize3 w e (cadr e))))))
  566.  
  567. (defun linearize3 (w e x)
  568.   (let (w1)
  569.        (if (and (memq w '($inf $minf $infinity)) (safe-zerop x))
  570.        (merror "Undefined form 0*inf:~%~M" e))
  571.        (setq w (mul2 (simplifya x t) w))
  572.        (cond ((or (atom w) (getl (caar w) '($outative $linear))) (setq w1 1))
  573.          ((eq (caar w) 'mtimes)
  574.           (setq w1 (cons '(mtimes) nil))
  575.           (do ((w2 (cdr w) (cdr w2)))
  576.           ((null w2) (setq w1 (nreverse w1)))
  577.           (if (or (atom (car w2))
  578.               (not (getl (caaar w2) '($outative $linear))))
  579.               (setq w1 (cons (car w2) w1)))))
  580.          (t (setq w1 w)))
  581.        (if (and (not (atom w1)) (or (among '$inf w1) (among '$minf w1)))
  582.        (infsimp w)
  583.        w)))
  584.  
  585. (setq opers (cons '$additive opers)
  586.       *opers-list (cons '($additive . additive) *opers-list))
  587.  
  588. (defun rem-opers-p (p)
  589.        (cond ((eq (caar opers-list) p)
  590.           (setq opers-list (cdr p)))
  591.          ((do ((l opers-list (cdr l)))
  592.           ((null l))
  593.           (if (eq (caar (cdr l)) p)
  594.               (return (rplacd l (cddr l))))))))
  595.  
  596. (defun additive (e z)
  597.   (cond ((get (caar e) '$outative)            ;Really a linearize!
  598.      (setq opers-list (copy-top-level opers-list ))
  599.      (rem-opers-p '$outative)
  600.      (linearize1 e z))
  601.     ((mplusp (cadr e))
  602.      (addn (mapcar #'(lambda (q)
  603.                  (oper-apply (list* (car e) q (cddr e)) z))
  604.                (cdr (cadr e)))
  605.            z))
  606.     (t (oper-apply e z))))
  607.  
  608. (setq opers (cons '$multiplicative opers)
  609.       *opers-list (cons '($multiplicative . multiplicative) *opers-list))
  610.  
  611. (defun multiplicative (e z)
  612.   (cond ((mtimesp (cadr e))
  613.      (simptimes
  614.       (cons '(mtimes)
  615.         (mapcar #'(lambda (q)
  616.                   (oper-apply (list* (car e) q (cddr e)) z))
  617.             (cdr (cadr e))))
  618.       1 z))
  619.     (t (oper-apply e z))))
  620.  
  621. (setq opers (cons '$outative opers)
  622.       *opers-list (cons '($outative . outative) *opers-list))
  623.  
  624. (defun outative (e z)
  625.        (setq e (cons (car e) (mapcar #'(lambda (q) (simpcheck q z)) (cdr e))))
  626.        (cond ((get (caar e) '$additive)
  627.           (setq opers-list (copy-top-level opers-list ))
  628.           (rem-opers-p '$additive)
  629.           (linearize1 e t))
  630.          ((linearconst e))
  631.          ((mtimesp (cadr e))
  632.           (let ((u (if (and (cddr e)
  633.                 (or (atom (caddr e))
  634.                     ($subvarp (caddr e))))
  635.                 (partition (cadr e) (caddr e) 1)
  636.                (rem-const (cadr e))))
  637.             (w))
  638.         (setq w (oper-apply (list* (car e)
  639.                        (simplifya (cdr u) t)
  640.                        (cddr e))
  641.                     t))
  642.         (linearize3 w e (car u))))
  643.          (t (oper-apply e t))))
  644.  
  645. (defprop %sum t $outative)
  646. (defprop %sum t opers)
  647. (defprop %integrate t $outative)
  648. (defprop %integrate t opers)
  649. (defprop %limit t $outative)
  650. (defprop %limit t opers)
  651.  
  652. (setq opers (cons '$evenfun opers)
  653.       *opers-list (cons '($evenfun . evenfun) *opers-list))
  654.  
  655. (setq opers (cons '$oddfun opers)
  656.       *opers-list (cons '($oddfun . oddfun) *opers-list))
  657.  
  658. (defmfun evenfun (e z)
  659.  (if (or (null (cdr e)) (cddr e))
  660.      (merror "Even function called with wrong number of arguments:~%~M" e))
  661.  (let ((arg (simpcheck (cadr e) z)))
  662.    (oper-apply (list (car e) (if (mminusp arg) (neg arg) arg)) t)))
  663.  
  664. (defmfun oddfun (e z)
  665.  (if (or (null (cdr e)) (cddr e))
  666.      (merror "Odd function called with wrong number of arguments:~%~M" e))
  667.  (let ((arg (simpcheck (cadr e) z)))
  668.    (if (mminusp arg) (neg (oper-apply (list (car e) (neg arg)) t))
  669.              (oper-apply (list (car e) arg) t))))
  670.  
  671. (setq opers (cons '$commutative opers)
  672.       *opers-list (cons '($commutative . commutative1) *opers-list))
  673.  
  674. (setq opers (cons '$symmetric opers)
  675.       *opers-list (cons '($symmetric . commutative1) *opers-list))
  676.  
  677. (defmfun commutative1 (e z)
  678.        (oper-apply (cons (car e)
  679.              (reverse
  680.               (sort (mapcar #'(lambda (q) (simpcheck q z))
  681.                     (cdr e))
  682.                 'great)))
  683.            t))
  684.  
  685. (setq opers (cons '$antisymmetric opers)
  686.       *opers-list (cons '($antisymmetric . antisym) *opers-list))
  687.  
  688. (declare-top (special sign))
  689.  
  690. (defmfun antisym (e z)
  691.  (let ((l (mapcar #'(lambda (q) (simpcheck q z)) (cdr e))))
  692.       (let (sign) (if (or (not (eq (caar e) 'mnctimes)) (freel l 'mnctimes))
  693.               (setq l (bbsort1 l)))
  694.           (cond ((equal l 0) 0)
  695.             ((prog1 (null sign) (setq e (oper-apply (cons (car e) l) t)))
  696.              e)
  697.             (t (neg e))))))
  698.  
  699. (defun bbsort1 (l)
  700.        (prog (sl sl1)
  701.          (if (or (null l) (null (cdr l))) (return l))
  702.          (setq sign nil sl (list nil (car l)))
  703.     loop (setq l (cdr l))
  704.          (if (null l) (return (nreverse (cdr sl))))
  705.          (setq sl1 sl)
  706.     loop1(cond ((null (cdr sl1)) (rplacd sl1 (cons (car l) nil)))
  707.            ((alike1 (car l) (cadr sl1)) (return 0))
  708.            ((great (car l) (cadr sl1)) (rplacd sl1 (cons (car l) (cdr sl1))))
  709.            (t (setq sign (not sign) sl1 (cdr sl1)) (go loop1)))
  710.          (go loop)))
  711.  
  712. (setq opers (cons '$nary opers)
  713.       *opers-list (cons '($nary . nary1) *opers-list))
  714.  
  715. (defmfun nary1 (e z)
  716.  (do ((l (cdr e) (cdr l)) (ans))
  717.      ((null l) (oper-apply (cons (car e) (nreverse ans)) z))
  718.      (setq ans (if (and (not (atom (car l))) (eq (caaar l) (caar e)))
  719.            (nconc (reverse (cdar l)) ans)
  720.            (cons (car l) ans)))))
  721.  
  722. (setq opers (cons '$lassociative opers)
  723.       *opers-list (cons '($lassociative . lassociative) *opers-list))
  724.  
  725. (defmfun lassociative (e z)
  726.   (let ((ans (cdr (oper-apply (cons (car e) (total-nary e)) z))))
  727.     (cond ((null (cddr ans)) (cons (car e) ans))
  728.       ((do ((newans (list (car e) (car ans) (cadr ans))
  729.             (list (car e) newans (car ans)))
  730.         (ans (cddr ans) (cdr ans)))
  731.            ((null ans) newans))))))
  732.  
  733. (setq opers (cons '$rassociative opers)
  734.       *opers-list (cons '($rassociative . rassociative) *opers-list))
  735.  
  736. ;(defmfun rassociative (e z)
  737. ;  (let ((ans (nreverse (cdr (oper-apply (cons (car e) (total-nary e)) z)))))
  738. ;    (cond ((null (cddr ans)) (cons (car e) ans))
  739. ;      ((do ((newans (list (car e) (cadr ans) (car ans))
  740. ;            (list (car e) (car ans) newans))
  741. ;        (ans (cddr ans) (cdr ans)))
  742. ;           ((null ans) newans))))))
  743. ;Update from F302 --gsb
  744. (defmfun rassociative (e z)
  745.   (let ((ans (cdr (oper-apply (cons (car e) (total-nary e)) z))))
  746.     (cond ((null (cddr ans)) (cons (car e) ans))
  747.       (t (setq ans (nreverse ans))
  748.          (do ((newans (list (car e) (cadr ans) (car ans))
  749.               (list (car e) (car ans) newans))
  750.           (ans (cddr ans) (cdr ans)))
  751.          ((null ans) newans))))))
  752.  
  753. (defmfun total-nary (e)
  754.   (do ((l (cdr e) (cdr l)) (ans))
  755.       ((null l) (nreverse ans))
  756.       (setq ans (if (and (not (atom (car l))) (eq (caaar l) (caar e)))
  757.             (nconc (reverse (total-nary (car l))) ans)
  758.             (cons (car l) ans)))))
  759.  
  760. #-NIL
  761. (setq opers (purcopy opers) *opers-list (purcopy *opers-list))
  762.  
  763. (defparameter $opproperties (cons '(mlist simp) (reverse opers)))
  764.  
  765.