src/lesfac.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 1981 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module lesfac)
  14
  15 (load-macsyma-macros rzmac ratmac)
  16
  17 (defun newsym2 (p e &aux (g (gensym)))
  18   (putprop g e 'disrep)
  19   (valput g (1- (valget (car genvar))))
  20   (setq genvar (cons g genvar))
  21   (setq varlist (cons e varlist))
  22   (putprop g p 'unhacked)
  23   g)
  24
  25 (defun getunhack (gen) (or (get gen 'unhacked) (pget gen)))
  26
  27 (defmacro getdis (x) `(get ,x 'disrep))
  28
  29 (defmacro cons1 (x) `(cons ,x 1))
  30
  31 (defun frpoly? (r) (equal 1 (cdr r)))
  32
  33 (defmacro setcall (&rest l)
  34   (setq l (cons 'setcall l))
  35   (sublis (list (cons 'fncall (cdr l))
  36                 (cons 'a (caddr l))
  37                 (cons 'b (cadddr l)))
  38           '(prog1 (car (setq a fncall)) (setq b (caddr a) a (cadr a)))))
  39
  40 (defun pquocof (p q)
  41   (let ((qq (testdivide p q)))
  42     (cond (qq (list q qq 1))
  43           ((list 1 p q)))))
  44
  45 (defun polyst (a)
  46   (cond ((pcoefp a) (list a))
  47         (t (cons (cons (car a) (cadr a)) (polyst (caddr a))))))
  48
  49 (defun cdinf (a b both)
  50   (cond ((or (pcoefp a) (pcoefp b)) (list 1 a b))
  51         (t (setq a (ncons (copy-tree a))
  52                  b (ncons (cond (both (copy-tree b))(t b))))
  53            (list (cd1 a b both) (car a) (car b)))))
  54
  55 (defun cd1 (a b both)
  56   (cond ((or (pcoefp (car a)) (pcoefp (car b))) 1)
  57         ((eq (caar a) (caar b))
  58          (ptimes (pexpt (pget (caar a)) ;CHECK FOR ALG. TRUNC.
  59                         (prog1 (cond (both (+ (cadar a) (cadar b))) (t (cadar a)))
  60                           (rplaca a (caddar a))
  61                           (cond (both (rplaca b (caddar b)))
  62                                 (t (setq b (cddar b))))))
  63                  (cd1 a b both)))
  64         ((pointergp (caar a) (caar b)) (cd1 (cddar a) b both))
  65         (t (cd1 a (cddar b) both))))
  66
  67 (defun lmake (p l)
  68   (cond ((pcoefp p) (cons p l))
  69         ((get (car p) 'unhacked)
  70          (lmake (caddr p) (cons (cons (car p) (cadr p)) l)))
  71         (t (setq l (lmake (caddr p) l))
  72            (rplaca l (list (car p) (cadr p) (car l))))))
  73
  74 (defun lmake2 (p l)
  75   (setq l (lmake p l))
  76   (mapc #'(lambda (x) (rplaca x (getunhack (car x))))
  77         (cdr l))
  78   (cond ((equal (car l) 1) (cdr l))
  79         (t (rplaca l (cons (car l) 1)))))
  80
  81
  82 (defun pmake (l)
  83   (cond ((null l) 1)
  84         ((= 0 (cdar l)) (pmake (cdr l)))
  85         ((numberp (caar l))          ;CLAUSE SHOULD BE ELIMINATED ASAP
  86          (ptimes (cexpt (caar l) (cdar l)) (pmake (cdr l))))
  87         (t (ptimes (list (caar l) (cdar l) 1) (pmake (cdr l))))))
  88
  89 (defun facmgcd (pl)            ;GCD OF POLY LIST FOR EZGCD WITH RATFAC
  90   (do ((l (cdr pl) (cdr l))
  91        (ans nil (cons (caddr gcd) ans))
  92        (gcd (car pl) (car gcd)))
  93       ((null l) (cons gcd (nreverse ans)))
  94     (setq gcd (fpgcdco gcd (car l)))
  95     (cond ((equal (car gcd) 1) (return (cons 1 pl)))
  96           ((null ans) (setq ans (list (cadr gcd))))
  97           ((not (equal (cadr gcd) 1))
  98            (do ((l2 ans (cdr l2))) ((null l2))
  99              (rplaca l2 (ptimes (cadr gcd) (car l2))))))))
 100
 101
 102 (defun fpgcdco (p q)
 103   (let ($ratfac gcdl)                   ;FACTORED PGCDCOFACTS
 104     (cond ((or (pcoefp p) (pcoefp q)) (pgcdcofacts p q))
 105           (t (list (ptimeschk
 106                     (setcall pgcdcofacts p q)
 107                     (car (setq p (lmake p nil)
 108                                q (lmake q nil)
 109                                gcdl (mapcar 'pmake (lgcd1 (cdr p) (cdr q)) ))))
 110                    (ptimeschk (car p) (cadr gcdl))
 111                    (ptimeschk (car q) (caddr gcdl)))))))
 112
 113 ;;      NOTE: ITEMS ON VARLIST ARE POS. NORMAL
 114 ;;      INTEGER COEF GCD=1 AND LEADCOEF. IS POS.
 115
 116 (defun lgcd1 (a b)
 117   (prog (ptlist g bj c t1 d1 d2)
 118      (setq ptlist (mapcar #'(lambda (ig) (declare (ignore ig)) b) a))
 119      (do ((a a (cdr a))
 120           (ptlist ptlist (cdr ptlist)))
 121          ((null a))
 122        (do ((ai (getunhack (caar a)))
 123             (b (car ptlist) (cdr b)))
 124            ((null b))
 125          (and (zerop (cdar b)) (go nextb))
 126          (setq d1 1 d2 1)
 127          (setq bj (getunhack (caar b)))
 128          (setq c (cond ((pirredp (caar a))
 129                         (cond ((pirredp (caar b)) 1)
 130                               (t (setcall pquocof bj ai))))
 131                        ((pirredp (caar b)) (setcall pquocof ai bj))
 132                        (t (setcall pgcdcofacts ai bj))))
 133          (cond ((equal c 1) (go nextb))
 134                ((equal ai 1) (go bloop)))
 135         aloop
 136          (cond ((setq t1 (testdivide ai c))
 137                 (setq ai t1 d1 (1+ d1))
 138                 (go aloop)))
 139         bloop
 140          (and (= d1 1)
 141               (not (equal bj 1))
 142               (do ((t1
 143                     (testdivide bj c)
 144                     (testdivide bj c)))
 145                   ((null t1))
 146                 (setq bj t1 d2 (1+ d2))))
 147          (setq g (cons (cons (makprodg c t)
 148                              (min (setq d1 (* d1 (cdar a)))
 149                                   (setq d2 (* d2 (cdar b)))))
 150                        g))
 151          (cond ((> d1 (cdar g))
 152                 (rplacd (last a)
 153                         (ncons (cons (caar g) (- d1 (cdar g)))))
 154                 (rplacd (last ptlist) (ncons (cdr b)))))
 155          (cond ((> d2 (cdar g))
 156                 (rplacd (last b)
 157                         (ncons (cons (caar g) (- d2 (cdar g)))))))
 158          (rplaca (car a) (makprodg ai t))
 159          (rplaca (car b) (makprodg bj t))
 160          (and (equal bj 1) (rplacd (car b) 0))
 161          (and (equal ai 1) (rplacd (car a) 0) (return nil))
 162         nextb))
 163      (return (list g a b))))
 164
 165 (defun makprodg (p sw)
 166   (cond ((pcoefp p) p)
 167         (t (car (makprod p sw)))))
 168
 169 (defun dopgcdcofacts (x y)
 170   (let (($gcd
 171          $gcd)( $ratfac nil)) (or (member $gcd *gcdl* :test #'eq) (setq $gcd '$ez))
 172          (pgcdcofacts x y)))
 173
 174 (defun facrplus (x y)
 175   (let ((a (car x))
 176         (b (cdr x))
 177         (c (car y))
 178         (d (cdr y)))
 179     (setq x (setcall dopgcdcofacts a c)
 180           y (setcall fpgcdco b d))
 181     (setq a (makprod
 182              (pplus (pflatten (ptimeschk a d))
 183                     (pflatten (ptimeschk b c))) nil))
 184     (setq b (ptimeschk b d))
 185     (cond ($algebraic (setq y (ptimeschk y b))
 186                       (setcall fpgcdco y a) ;for unexpected gcd
 187                       (cons (ptimes x a) y))
 188           (t (setq c (setcall cdinf y b nil))
 189              (setcall fpgcdco y a)
 190              (cons (ptimes x a) (ptimeschk y (ptimeschk c b)))))))
 191
 192 (defun mfacpplus (l)
 193   (let (($gcd (or $gcd '$ez))
 194         ($ratfac nil)
 195         (g nil))
 196     (setq g (oldcontent2 (sort (copy-list l) 'contodr) 0))
 197     (cond ((pzerop g) g)
 198           ((do ((a (pflatten (pquotient (car l) g))
 199                    (pplus a (pflatten (pquotient (car ll) g))))
 200                 (ll (cdr l) (cdr ll)))
 201                ((null ll) (ptimes g (makprod a nil))))))))
 202
 203 (defun  facrtimes (x y gcdsw)
 204   (cond ((not gcdsw)
 205          (cons (ptimes (car x) (car y)) (ptimeschk (cdr x) (cdr y))))
 206         (t (let ((g (cdinf (car x) (car y) t))
 207                  (h (cdinf (cdr x) (cdr y) t)))
 208              (setq x (fpgcdco (cadr g) (caddr h)))
 209              (setq y (fpgcdco (caddr g) (cadr h)))
 210              (cons (ptimes (car g) (ptimes (cadr x) (cadr y)))
 211                    (ptimeschk (car h) (ptimeschk (caddr x) (caddr y))))))))
 212
 213 (defun pfacprod (poly)                  ;FOR RAT3D
 214   (cond ((pcoefp poly) (cfactor poly))
 215         (t (nconc (pfacprod (caddr poly))
 216                   (list (pget (car poly)) (cadr poly))))))
 217
 218 (defun fpcontent (poly)
 219   (let (($ratfac
 220          nil))                          ;algebraic uses
 221     (setq poly (oldcontent poly))       ;rattimes?
 222     (let ((a (lowdeg (cdadr poly))))    ;main var. content
 223       (cond ((> a 0) (setq a (list (caadr poly) a 1))
 224              (setq poly
 225                    (list (ptimes (car poly) a)
 226                          (pquotient (cadr poly) a))))))
 227     (if (pminusp (cadr poly))
 228         (list (pminus (car poly)) (pminus (cadr poly)))
 229         poly)))
 230
 231 ;; LOWDEG written to compute the lowest degree of a polynomial. - RZ
 232
 233 (defmfun lowdeg (p)
 234   (do ((l p (cddr l)))
 235       ((null (cddr l)) (car l))))
 236
 237 (defun makprod (poly contswitch)
 238   (cond ((pureprod poly) poly)
 239         ((null (cdddr poly))
 240          (ptimes (list (car poly) (cadr poly) 1)
 241                  (makprod (caddr poly) contswitch)))
 242         (contswitch (makprod1 poly))
 243         (t (setq poly (fpcontent poly))
 244            (ptimes (makprod (car poly) contswitch) (makprod1 (cadr poly))))))
 245
 246 (defun makprod1 (poly)
 247   (do ((v varlist (cdr v))
 248        (g genvar (cdr g))
 249        (p (pdis poly)))
 250       ((null v) (maksymp poly))
 251     (and (alike1 p (car v)) (return (pget (car g))))))
 252
 253 (defun maksym (p)
 254   (newsym2 p (pdis p)))
 255
 256 (defun maksymp (p)
 257   (cond ((atom p) p)
 258         (t (pget (maksym p)))))
 259
 260 (defun pflatten (h)
 261   (prog (m)
 262      (setq m (listovars h))
 263      checkmore
 264      (cond ((null m) (return h))
 265            ((not (let ((p (getunhack (car m))))
 266                    (or (null p) (eq (car m) (car p)))))
 267             (go redo))
 268            (t (setq m (cdr m)) (go checkmore)))
 269      redo (return (let ($ratfac) (pflat1 h)))))
 270
 271 (defun pflat1 (p)
 272   (cond ((pcoefp p) p)
 273         ((null (cdddr p))
 274          (ptimes (pexpt (getunhack (car p)) (cadr p)) (pflat1 (caddr p))))
 275         (t (do ((val (getunhack (car p)))
 276                 (ld (cadr p) (car a))
 277                 (a (cdddr p) (cddr a))
 278                 (ans (pflat1 (caddr p))))
 279                ((null a) (ptimes ans (pexpt val ld)))
 280              (setq ans
 281                    (pplus (ptimes ans
 282                                   (pexpt val (- ld (car a))))
 283                           (pflat1 (cadr a))))))))
 284
 285 (defun pirredp (x)
 286   (and (setq x (getdis x))
 287        (or (atom x) (member 'irreducible (cdar x) :test #'eq))))
 288
 289 (defun knownfactors (d)
 290   (prog (h)
 291      (cond ((pcoefp d) (return d)))
 292      (setq h (getdis (car d)))
 293      (return (cond ((or (atom h) (not (eq (caar h) 'mtimes)))
 294                     (ptimes (knownfactors (caddr d))
 295                             (list (car d) (cadr d) 1)))
 296                    (t (setq h (getunhack (car d)))
 297                       (ptimes (knownfactors (caddr d))
 298                               (pexpt (knownfactors h) (cadr d))))))))