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