src/ufact.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12 (macsyma-module ufact)
  13
  14 (load-macsyma-macros ratmac rzmac)
  15
  16 ;; Dense Polynomial Representation
  17
  18 (defun dprep (p)
  19   (do ((n (car p))
  20        (e (car p) (f1- e))
  21        (l))
  22       ((< e 0) (cons n (nreverse l)))
  23     (cond ((equal e (car p))
  24            (push (cadr p) l)
  25            (setq p (cddr p)))
  26           (t (push 0 l)))))
  27
  28 (defun dpdisrep (l)
  29   (cond ((zerop (car l)) (cadr l))
  30         ((do ((l (nreverse (cdr l)) (cdr l))
  31               (n 0 (f1+ n))
  32               (ll))
  33              ((null l) ll)
  34            (or (= (car l) 0)
  35                (setq ll (cons n (cons (car l) ll))))))))
  36
  37 ;; not currently called
  38 ;;(DEFUN PGCDU* (P Q)
  39 ;;       (COND ((OR (PCOEFP P) (PCOEFP Q))  1)
  40 ;;           ((NULL MODULUS)
  41 ;;            (merror  "Illegal CALL TO PGCDU"))
  42 ;;           ((> (CADR P) (CADR Q))
  43 ;;            (PSIMP (CAR P) (DPDISREP (DPGCD (DPREP (CDR P)) (DPREP (CDR Q))))))
  44 ;;           ((PSIMP (CAR P) (DPDISREP (DPGCD (DPREP (CDR Q)) (DPREP (CDR P))))))))
  45 ;;
  46 ;;(DEFUN PMODSQFRU (P)
  47 ;;       (DO ((DPL (DPSQFR (DPREP (CDR P))) (CDR DPL))
  48 ;;          (PL NIL (CONS (PSIMP (CAR P) (DPDISREP (CDAR DPL))) (CONS (CAAR DPL) PL))))
  49 ;;         ((NULL DPL) PL)))
  50 \f
  51 (defun dpgcd (p q)
  52   (if (< (car p) (car q)) (rotatef p q))
  53   (do ((p (copy-list p) q)
  54        (q (copy-list q) (dpremquo p q nil)))
  55       ((= (car q) 0)
  56        (if (= (cadr q) 0) p '(0 1)))))
  57
  58 (defun dpdif (p q)
  59   (cond ((> (car p) (car q))
  60          (do ((i (car p) (f1- i))
  61               (pl (cdr p) (cdr pl))
  62               (l nil (cons (car pl) l)))
  63              ((= i (car q)) (dpdif1 pl (cdr q) l)) ))
  64         ((< (car p) (car q))
  65          (do ((i (car q) (f1- i))
  66               (ql (cdr q) (cdr ql))
  67               (l nil (cons (cminus (car ql)) l)))
  68              ((= i (car p)) (dpdif1 (cdr p) ql l))))
  69         (t (dpdif1 (cdr p) (cdr q) nil))))
  70
  71 (defun dpdif1 (p1 q1 l)
  72   (do ((pl p1 (cdr pl))
  73        (ql q1 (cdr ql))
  74        (ll l (cons (cdifference (car pl) (car ql)) ll)))
  75       ((null pl) (dpsimp (nreverse ll)))))
  76
  77 (defun dpsimp (pl) (setq pl (ufact-strip-zeroes pl))
  78        (cond ((null pl) '(0 0))
  79              (t (cons (f1- (length pl)) pl))))
  80
  81 (defun dpderiv (p)
  82   (cond ((= 0 (car p)) '(0 0))
  83         (t (do ((l (cdr p) (cdr l))
  84                 (i (car p) (f1- i))
  85                 (dp nil (cons (ctimes i (car l)) dp)))
  86                ((= i 0) (cons (f1- (car p)) (nreverse dp)))))))
  87
  88 (defun dpsqfr (q)                       ;ASSUMES MOD > DEGREE
  89   (do ((c q (dpmodquo c p))
  90        (d (dpderiv q) (dpmodquo d p))
  91        (i 0 (f1+ i))
  92        (p)
  93        (pl))
  94       ((= 0 (car c)) pl)
  95     (cond (p (setq d (dpdif d (dpderiv c))
  96                    p (dpgcd c d))
  97              (and (> (car p) 0)
  98                   (setq pl (cons (cons i p) pl))))
  99           (t (setq p (dpgcd c d))
 100              (cond ((= (car p) 0) (return (ncons (cons 1 c)))))))))
 101
 102
 103 \f
 104 (defun dpmodrem (p q)
 105   (cond ((< (car p) (car q)) p)
 106         ((= (car q) 0) '(0 0))
 107         ((dpremquo (copy-list p) (copy-list q) nil))))
 108
 109 (defun dpmodquo (p q)
 110   (cond ((< (car p) (car q)) '(0 0))
 111         ((= (car q) 0)
 112          (cond ((equal (cadr q) 1) p)
 113                (t (cons (car p)
 114                         (mapcar #'(lambda (c) (cquotient c (cadr q))) (cdr p))
 115                         ))))
 116         ((dpremquo (copy-list p) (copy-list q) t))))
 117
 118 ;; If FLAG is T, return quotient.  Otherwise return remainder.
 119
 120 (defun dpremquo (p q flag)
 121   (prog (lp lq l alpha)
 122      (cond ((= (cadr q) 1)
 123             (setq alpha 1))
 124            (t (setq alpha (crecip (cadr q)))
 125               (do ((l (cddr q) (cdr l)))
 126                   ((null l)
 127                    (rplaca (cdr q) 1))
 128                 (rplaca l (ctimes (car l) alpha)))))
 129      a    (and flag (setq l (cons (ctimes (cadr p) alpha) l)))
 130      (setq lp (cddr p) lq (cddr q))
 131      b    (rplaca lp (cdifference (car lp) (ctimes (car lq) (cadr p))))
 132      (cond ((null (setq lq (cdr lq)))
 133             (do ((e (f1- (car p)) (f1- e))
 134                  (pp (cddr p) (cdr pp)))
 135                 ((null pp) (setq p '(0 0)))
 136               (cond ((signp e (car pp))
 137                      (and flag (not (< e (car q)))
 138                           (setq l (cons 0 l))))
 139                     ((return (setq p (cons e pp))))))
 140             (cond ((< (car p) (car q))
 141                    (return (cond (flag (dpsimp (nreverse l)));GET EXP?
 142                                  (p))))
 143                   ((go a))))
 144            (t (setq lp (cdr lp))
 145               (go b)))))
 146
 147 (defun ufact-strip-zeroes (l)
 148   (do ((l l (cdr l)))
 149       ((null (pzerop (car l))) l)))
 150
 151 (defun cpres1 (a b)
 152   (prog (res (v 0) a3) (declare (fixnum v))
 153         (setq  a (dprep a) b (dprep b))
 154         (setq res 1)
 155         again (setq a3 (dpmodrem a b))
 156         (setq v (boole boole-xor v (logand 1 (car a) (car b) )))
 157         (setq res (ctimes res (cexpt (cadr b)
 158                                      (f- (car a) (car a3)))))
 159         (cond ((= 0 (car a3))
 160                (setq res (ctimes res (cexpt (cadr a3) (car b))))
 161                (return (cond ((oddp v) (cminus res))
 162                              (t res))) ))
 163         (setq a b)
 164         (setq b a3)
 165         (go again) ))