src/sumcon.lisp

   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
  13 (macsyma-module sumcon)
  14
  15 (declare-top (special $genindex $niceindicespref $sumexpand))
  16
  17 (defmfun $sumcontract (e)              ; e is assumed to be simplified
  18   (cond ((atom e) e)
  19         ((eq (caar e) 'mplus)
  20          (do ((x (cdr e) (cdr x)) (sums) (notsums) (car-x))
  21              ((null x) (cond ((null sums)
  22                               (subst0 (cons '(mplus)
  23                                             (nreverse notsums))
  24                                       e))
  25                              (t (setq sums (sumcontract1 sums))
  26                                 (addn (cons sums notsums) t))))
  27            (setq car-x (car x))
  28            (cond ((atom car-x)
  29                   (setq notsums (cons car-x notsums)))
  30                  ((eq (caar car-x) '%sum)
  31                   (setq sums (cons (cons ($sumcontract (cadr car-x))
  32                                          (cddr car-x))
  33                                    sums)))
  34                  (t (setq notsums (cons car-x notsums))))))
  35         (t (recur-apply #'$sumcontract e))))
  36
  37 (defmfun $intosum (e)                  ; e is assumed to be simplified
  38   (let (($sumexpand t))
  39     (cond ((atom e) e)
  40           ((eq (caar e) 'mtimes)        ;puts outside product inside
  41            (do ((x (cdr e) (cdr x)) (sum) (notsum))
  42                ((null x) (cond ((null sum)
  43                                 (subst0 (cons '(mtimes)
  44                                               (nreverse notsum))
  45                                         e))
  46                                (t (simpsum
  47                                    (let ((new-index
  48                                           (if (free (cons nil notsum) (caddr sum))
  49                                               (caddr sum)
  50                                               (get-free-index (cons nil (cons sum notsum))))))
  51                                      (setq sum (subst new-index (caddr sum) sum))
  52                                      (rplaca (cdr sum) (muln (cons (cadr sum) notsum) t))
  53                                      (rplacd (car sum) nil)
  54                                      sum)
  55                                    1 t))))
  56              (cond ((atom (car x))
  57                     (setq notsum (cons (car x) notsum)))
  58                    ((eq (caaar x) '%sum)
  59                     (setq sum (if (null sum)
  60                                   (copy-tree (car x))
  61                                   (muln (list sum (car x)) t))))
  62                    (t (setq notsum (cons ($sumcontract (car x))
  63                                          notsum))))))
  64           (t (recur-apply #'$intosum e)))))
  65
  66 (defun sumcontract1 (sums)
  67   (addn (sumcontract2 nil sums) t))
  68
  69 (defun sumcontract2 (result left)
  70   (if (null left)
  71       result
  72       (let ((x (sumcombine1 (car left) (cdr left))))
  73         (sumcontract2 (append (car x) result) (cdr x)))))
  74
  75 (defun sumcombine1 (pattern llist)
  76   (do ((sum pattern) (non-sums nil)
  77        (un-matched-sums nil) (try-this-one)
  78        (llist llist (cdr llist)))
  79       ((null llist) (cons (cons (simplify (cons '(%sum) sum))
  80                                 non-sums)
  81                           un-matched-sums))
  82     (setq try-this-one (car llist))
  83     (cond ((and (numberp (sub* (caddr sum) (caddr try-this-one)))
  84                 (numberp (sub* (cadddr sum) (cadddr try-this-one))))
  85            (let ((x (sumcombine2 try-this-one sum)))
  86              (setq sum (cdar x)
  87                    non-sums (cons (cdr x) non-sums))))
  88           (t (setq un-matched-sums (cons try-this-one un-matched-sums))))))
  89
  90 (defun sumcombine2 (sum1 sum2)
  91   (let* ((e1 (car sum1))
  92          (e2 (car sum2))
  93          (i1 (cadr sum1))
  94          (i2 (cadr sum2))
  95          (l1 (caddr sum1))
  96          (l2 (caddr sum2))
  97          (h1 (cadddr sum1))
  98          (h2 (cadddr sum2))
  99          (newl (simplify `(($max) ,l1 ,l2)))
 100          (newh (simplify `(($min) ,h1 ,h2)))
 101          (newi (cond ((eq i1 i2) i1)
 102                      ((free e1 i2) i2)
 103                      ((free e2 i1) i1)
 104                      (t (get-free-index (list nil i1 i2 e1 e2   l1 l2 h1 h2)))))
 105          (extracted nil)
 106          (new-sum nil))
 107     (setq e1 (subst newi i1 e1))
 108     (setq e2 (subst newi i2 e2))
 109     (setq new-sum (list '(%sum) (add2 e1 e2) newi newl newh))
 110     (setq extracted
 111           (addn
 112            (mapcar #'dosum
 113                    (list e1 e1 e2 e2)
 114                    (list newi newi newi newi)
 115                    (list l1 (add2 newh 1) l2 (add2 newh 1))
 116                    (list (sub* newl 1) h1 (sub* newl 1) h2)
 117                    '(t t t t))
 118            t))
 119     (cons new-sum extracted)))
 120
 121 (defmvar $niceindicespref '((mlist simp) $i $j $k $l $m $n))
 122
 123 (putprop '$niceindicespref 'assign-nonempty-list 'assign)
 124
 125 (defun assign-nonempty-list (n v)
 126   (declare (ignore n))
 127   (unless (and ($listp v) (not ($emptyp v)))
 128     (merror
 129       (intl:gettext "niceindicespref: value must be a nonempty list; found: ~:M")
 130       v)))
 131
 132 (defun get-free-index (llist &optional i)
 133   (or (do ((try-list (cdr $niceindicespref) (cdr try-list)))
 134           ((null try-list))
 135         (if (or (free llist (car try-list))
 136                 (eq i (car try-list)))
 137             (return (car try-list))))
 138       (do ((n 0 (1+ n)) (try))
 139           (nil)
 140         (setq try (intern (format nil "~a~d" (cadr $niceindicespref) n)))
 141         (if (free llist try) (return try)))))
 142
 143 (defmfun $bashindices (e)              ; e is assumed to be simplified
 144   (if (atom e)
 145     e
 146     (let (($genindex '$j)
 147           (e (recur-apply #'$bashindices e)))
 148       (cond ((atom e) e)
 149             ((member (caar e) '(%sum %product) :test #'eq)
 150              (sumconsimp (subst (gensumindex) (caddr e) e)))
 151             (t e)))))
 152
 153 (defmfun $niceindices (e)
 154   (if (atom e)
 155       e
 156       (let ((e (recur-apply #'$niceindices e)))
 157         (cond ((atom e) e)
 158               ((member (caar e) '(%sum %product) :test #'eq)
 159                (sumconsimp (subst (get-free-index e (caddr e)) (caddr e) e)))
 160               (t e)))))
 161
 162 (defun sumconsimp (e)
 163   (if (and (not (atom e)) (member (caar e) '(%sum %product) :test #'eq))
 164       (list* (car e) (sumconsimp (cadr e)) (cddr e))
 165       (resimplify e)))