share/affine/ndotsimp.lisp

   1 ;;; -*- Package: CL-MAXIMA; Mode: LISP; Syntax: Common-lisp -*-
   2 (in-package :maxima)
   3
   4 ;(defun $check_associative (a b c &aux tem tem1 tem2)
   5 ;;  (show (dotsimp (n. a b)))
   6 ;;  (show (dotsimp (n. b c)))
   7 ;;  (show (n. a (dotsimp (n. b c))))
   8 ;;  (show (setq hee(n.  (dotsimp (n. a b)) c)))
   9 ;  (setq tem1 (new-dotsimp(n. (new-dotsimp (n. a b)) c)))
  10 ;  (setq tem2 (new-dotsimp(n. a (new-dotsimp (n. b c)))))
  11 ;;  (show tem1 tem2)
  12 ;  (setq tem (n- tem1 tem2))
  13 ;  (cond ((pzerop tem) tem)
  14 ;       (t (header-poly tem))))
  15
  16 (defun $check_associative (a b c &aux tem tem1 tem2)
  17   (setq tem1 (n. (new-dotsimp (n. a b)) c))
  18   (setq tem2 (n. a (new-dotsimp (n. b c))))
  19   (setq tem (n- tem1 tem2))
  20   (setq tem (new-dotsimp tem))
  21   (cond ((pzerop tem) tem)
  22         (t (header-poly tem))))
  23
  24 (defun ncdot-list (lis)
  25   (cond ((null lis) 1)
  26         ((null (cdr lis))
  27          (car lis))
  28         (t (cons '(mnctimes) lis))))
  29
  30 (defun dot-subword (small big &aux tem leng)
  31   "(dot-subword #$x$ #$z.u.x.y.z.z.y$) ==> (list #$z.u$ #$y.z.z.y$))
  32   (dot-subword #$x.y$ #$z.u.x.y.z.z.y$) ==> (list #$z.u$ #$z.z.y$))
  33   (dot-subword #$z.u.x.y$ #$z.u.x.y.z.z.y$) ==> (list #$1$ #$z.z.y$))"
  34   (cond ((atom big)
  35          (cond ((eql big small) (list 1  1))))
  36         ((atom small)
  37          (cond ((setq tem (member small big :test #'eq))
  38                 (list (ncdot-list (subseq (cdr big) 0 (- (length big) (length tem) 1)))
  39                       (ncdot-list (cdr tem))))))
  40         (t (setq leng (-  (length small) 1))
  41            (loop for v on (cdr big)
  42                  when ;;first part of v is equal to small
  43                  (loop
  44                    initially (cond ((> leng (length v)) (return nil)))
  45                    for vv in v
  46                    for w in (cdr small)
  47                    when (not (equal w vv))
  48                    do (return nil)
  49                    finally (return t))
  50                  do (return (list  (ncdot-list
  51                                      (subseq (cdr big) 0 (- (length big) (length v) 1)))
  52                                    (ncdot-list (nthcdr (length (cdr small)) v))))))))
  53
  54 (defun split-numerator (reduced rest)
  55   "splits into two polynomials, the first one needs no replacement and the second has its leading term needing
  56    replacement"
  57   (cond ((pzerop rest) (values reduced 0))
  58         ((poly-scalarp rest) (values (n+ rest reduced) 0))
  59         (($must_replacep (get (p-var rest) 'disrep))
  60          (values reduced rest))
  61         (t (split-numerator  (n+ (subseq rest 0 3) reduced) (or (fifth rest) 0)))))
  62
  63 (defun new-dotsimp (ratl-fun &aux mon repl num den)
  64   (format t "~%Beginning to simplify:")
  65   (sh ratl-fun)
  66   (loop
  67      with expr = ratl-fun with answer = 0
  68      when (pzerop expr) do (return answer)
  69      do   (setq-num-den num den expr)
  70        (cond ((poly-scalarp num)(setq answer (n+ answer expr))
  71               (format t "~%Final answer:")
  72               (sh answer)
  73               (return answer))
  74              (($must_replacep (setq mon (get (p-var num) 'disrep)))
  75               (cond (*verbose-check-overlaps*
  76                      (format t "~%Simplifying the worst monomial ")
  77                      (dot-show mon)))
  78               (setq repl (simp-once-monomial mon))
  79               (setq repl (n* (p-cof num) repl))
  80               (cond ((fifth num)
  81                      (setq num (n+ repl (fifth num))))
  82                     (t (setq num repl)))
  83               (setq expr (nred num den)))
  84              (t (multiple-value-bind
  85                       (reduced rest)
  86                     (split-numerator 0 num)
  87                   (setq answer (n+ (nred reduced den) answer))
  88                   (setq expr (nred rest den))
  89                   (cond (*verbose-check-overlaps*
  90                          (format t "~%Simplifying the worst monomial ")
  91                          (dot-show mon) (format t " adding  to the answer" )))
  92                   )))
  93                                         ;      (format t "~%Expr:") (sh expr)
  94                                         ;      (format t "~%Reduced part:") (sh answer)
  95        ))
  96
  97 ;;;old reliable
  98 ;(defun new-dotsimp (ratl-fun &aux mon repl num den)
  99 ;  (progn
 100 ;    (loop
 101 ;      with expr = ratl-fun with answer = 0
 102 ;      when (pzerop expr) do (return answer)
 103 ;      do   (setq-num-den num den expr)
 104 ;      (cond ((poly-scalarp num)(setq answer (n+ answer expr))
 105 ;                              (return answer))
 106 ;           (($must_replacep (setq mon (get (p-var num) 'disrep)))
 107 ;            (format t "~%Simplifying the worst monomial ") (dot-show mon)
 108 ;            (setq repl (simp-once-monomial mon))
 109 ;            (setq repl (n* (p-cof num) repl))
 110 ;            (cond ((fifth num)
 111 ;                   (setq num (n+ repl (fifth num))))
 112 ;                  (t (setq num repl)))
 113 ;            (setq expr (nred num den)))
 114 ;           (t (setq answer (n+ (nred (firstn 3 num) den) answer))
 115 ;              (format t "~%Simplifying the worst monomial ")
 116 ;              (dot-show mon) (format t " adding to the answer")
 117 ;              (cond ((fifth num)
 118 ;                     (setq expr (nred (fifth num) den)))
 119 ;                    (t (setq expr 0)))))
 120 ;      )))
 121
 122 (defun force-poly (repl)
 123   (cond ((numberp repl) repl)
 124         (t (cdr repl))))
 125
 126 (defvar $dot_eps nil)
 127
 128 (defun simp-once-monomial (monom &aux tem)
 129   (cond ((atom monom)
 130          (loop for (mon repl) on (cdr $dot_simplifications) by #'cddr
 131                when (eql mon monom)
 132                do (return (force-poly repl))))
 133         ((and $dot_eps
 134              (setq tem (member '$eps monom :test 'eq))
 135              (member'$eps (cdr tem) :test 'eq))
 136          0)
 137         (t (loop for (mon repl) on (cdr $dot_simplifications) by #'cddr
 138                 when (setq tem (dot-subword mon monom))
 139                 do
 140                 (return (n. (n.  (first tem) (force-poly repl)) (second tem)))
 141                 finally (return (st-rat monom))))))
 142
 143 (defun $dotsimp (expr)
 144   (declare (special $new_fast_dotsimp))
 145   (cond (($listp expr) (cons '(mlist) (mapcar '$dotsimp (cdr expr))))
 146         (t
 147          (cond ((or (rational-functionp expr)(affine-polynomialp expr)) nil)
 148                (t (setq expr (new-rat expr))))
 149          (header-poly (new-dotsimp expr)))))
 150
 151
 152
 153
 154 (defun $dot_factor (form variables &optional (slot 1))
 155 "Form is a non commutative polynomial, and variables a list of non commutative variables.
 156 Form is split into a list of forms of length variables, so that
 157 variables.dot_factor(form) = form "
 158
 159   (check-arg variables $listp "macsyma list")
 160   (let* ((fo (st-rat form))
 161          (vari (list-variables fo))
 162          (result (make-list (length variables) :initial-element 0))
 163          pos vv)
 164     (loop for v in vari
 165           do (setq vv (get v 'disrep))
 166           (cond ((not (symbolp vv))
 167                  (setq pos (position (nth slot vv) variables))
 168                  (cond (pos (setf (nth pos result)
 169                                   (n+ (nth pos result) (n* (pcoeff fo v) (st-rat
 170                                                            (meval* (cons '(mnctimes)
 171                                                                      (let ((mon (copy-list (cdr vv))))
 172                                                                         (setf (nth (1- slot) mon) 1) mon))))))))
 173                        (t (fsignal "does not have a factor in slot"))))
 174                 ((setq pos (position vv variables))
 175                  (cond (pos (setf (nth pos result)
 176                                   (n+ (nth pos result) (pcoeff fo (list v 1 1)))))))))
 177     (cons '(mlist) (mapcar 'new-disrep (cdr result)))))