src/algfac.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
  13 (macsyma-module algfac)
  14
  15 ;; this is the alg factor package
  16
  17 ;;; Toplevel functions are: CPBGZASS CPTOM
  18
  19 (load-macsyma-macros ratmac)
  20
  21 (declare-top (special tra* trl* *xn var intbs* plim many* split* alc ind p l))
  22
  23 (defun ziredup (p)
  24   (let ((modulus nil) (*alpha* nil) (minpoly* nil) (algfac* nil)
  25         (gauss nil) (tellratlist nil) (many* nil)
  26         (mm* 1)
  27         ($gcd '$ez))
  28     (null (cddr(pfactor p)))))
  29
  30 (defun intbasehk (p)
  31   (prog (modulus)
  32      (setq modulus plim)
  33      (setq p (pctimes intbs* p))
  34      (setq modulus nil)
  35      (return (car (ratreduce p intbs*)))))
  36
  37 (defun findibase (p)
  38   (prog (mainvar)
  39      (setq mainvar (car p))
  40      (setq p (redresult p (pderivative p mainvar)))
  41      (setq p (cfactorw p))
  42      (setq mainvar 1)
  43      loop (when (null p) (return mainvar))
  44      (setq mainvar (* mainvar (expt (car p) (quotient (cadr p) 2))))
  45      (setq p (cddr p))
  46      (go loop)))
  47
  48
  49 (defun cpbgzass (qlist v m)
  50   (prog (f y vj factors u w lc j p2 fnj fnq oldfac)
  51      (cond ((equal m 1)
  52             (return (list v)))
  53            ((equal m (cadr v))
  54             (return (let ((var (list var 1 1)))
  55                       (gfsplit v)))))
  56      (setq f (pmod v))
  57      (setq lc (caddr f))
  58      (setq f (monize f))
  59      (setq p2 1
  60            qlist (cdr (nreverse qlist)))
  61      (setq oldfac (list nil f))
  62      nextq(setq v (car qlist))
  63      (setq qlist (cdr qlist))
  64      (setq j (findses v f))
  65      (setq oldfac (nconc oldfac fnq))
  66      (setq fnq nil)
  67      incrj(setq factors (nconc oldfac fnj))
  68      (setq fnj nil)
  69      (setq vj (pplus v (car j))
  70            j (cdr j))
  71      tag2 (setq u (cadr factors))
  72      (setq w (pgcdu vj u))
  73      (when (or (numberp w) (and *alpha* (alg w))(= (cadr u) (cadr w)))
  74        (go nextfac))
  75      (setq y (car (pmodquo u w)))
  76      (setq fnq (cons w fnq))
  77      (setq fnj (cons y fnj))
  78      (incf p2)
  79      (rplacd factors (cddr factors))
  80      (if (equal p2 m)
  81          (go out)
  82          (go tag1))
  83      nextfac
  84      (setq factors (cdr factors))
  85      tag1 (cond ((cdr factors)
  86                  (go tag2))
  87                 (j (go incrj))
  88                 (qlist (go nextq)))
  89      out  (setq fnq (nconc fnq fnj (cdr oldfac)))
  90      (return (cons (ptimes lc (car fnq)) (cdr fnq)))))
  91
  92
  93 ;; The function PMONZ used to be defined here.  It is also defined in
  94 ;; RAT;RAT3A and BMT claims the definitions are equivalent.
  95
  96 (defun findses (g f)
  97   (prog (var tra* trl*)
  98      (setq g (zassg (cdr g) (cdr f) (car g)))
  99      (setq var (list (car f) 1 1))
 100      (setq f (gfsplit g))
 101      (return (mapcar #'(lambda (a) (car (last a))) f))))
 102
 103 (defun coefvec (p n vec)
 104   (prog nil
 105    loop (when (zerop n) (return vec))
 106    (decf n)
 107    (push (ptterm p n) vec)
 108    (go loop)))
 109
 110 (defun zassg (g f var)
 111   (prog (i mat gn ans n)
 112      (setq n (car f))
 113      (setq gn g)
 114      (setq i 1
 115            mat (list (coefvec '(0 1) n (list 1))))
 116      (go on)
 117      loop (incf i)
 118      (setq gn (pgcd1 (ptimes1 gn g) f))
 119      on   (setq ans (lindep mat (coefvec gn n (list (list var i 1)))))
 120      (cond (ans (return ans)))
 121      (go loop)))
 122
 123 (defun divl (j a)
 124   (mapcar #'(lambda (l) (car (pmodquo l a))) j))
 125
 126 ;; (DEFUN PADDROWS (A B) (MAPCAR (FUNCTION PPLUS) A B))
 127
 128 (defun pdifrows (a b)
 129   (mapcar #'pdifference a b))
 130
 131 (defun ptimesrow (var row)
 132   (mapcar #'(lambda (a) (ptimes var a)) row))
 133
 134 (defun ddiv (j)
 135   (prog (a b)
 136      (setq b j)
 137      ag   (setq a (car b))
 138      (cond ((zerop a)
 139             (setq b (cdr b))
 140             (go ag)))
 141      (return (divl j a))))
 142
 143 (defun lindep (mat vec)
 144   (prog (e d m row rowd vecd)
 145      (setq m mat)
 146      (cond ((equal 0. (car vec)) (setq vec (cdr vec)))
 147            (t (setq vec (pdifrows (cdr vec) (ptimesrow (car vec) (cdar mat))))))
 148      loop (cond ((null (cdr m))
 149                  (cond ((zerolp (cdr (reverse vec)))
 150                         (return (car (last vec))))
 151                        (t (rplacd m (cons (ddiv vec) (cdr m)))
 152                           (return nil)))))
 153      (setq row (cadr m))
 154      (setq rowd row vecd vec)
 155      loop1(setq d (car rowd))
 156      (setq e (car vecd))
 157      (cond ((equal 0 e)
 158             (cond ((equal 0 d)
 159                    (setq vecd (cdr vecd) rowd (cdr rowd))
 160                    (go loop1))
 161                   (t (setq vec (cdr vec)) (setq m (cdr m)) (go loop))))
 162            ((equal 0 d)
 163             (rplacd m
 164                     (cons (divl vec e) (mapcar (function cdr) (cdr m))))
 165             (return nil)))
 166      (setq vec (pdifrows (cdr vec) (ptimesrow e (cdr row))))
 167      (setq m (cdr m))
 168      (go loop)))
 169
 170 (defun gfsplit (f)
 171   (prog (tr fl (n 0) ans tra* (i 0) nfl)
 172      (setq fl (list f) n (cadr f))
 173      loop (cond ((null fl)
 174                  (cond ((null nfl)
 175                         (cond ((= n (length ans))
 176                                (setq trl* nil)
 177                                (return ans))
 178                               (t (merror (intl:gettext "GFSPLIT: unknown error.")))))
 179                        (t
 180                         (setq fl nfl
 181                               nfl nil
 182                               i (1+ i))))))
 183      (setq f (car fl)
 184            fl (cdr fl))
 185      (cond ((> i mm*)
 186             (merror (intl:gettext "GFSPLIT: unknown error."))))
 187      (setq tr (tracemod0 f i))
 188      (cond ((or (pcoefp tr) (and algfac* (alg tr)))
 189             (setq nfl (cons f nfl))
 190             (go loop)))
 191      (setq f (cpbg0 tr f))
 192      (setq ans (nconc ans (car f)))
 193      (when (null (cdr f)) (go loop))
 194      (setq nfl (nconc nfl (cdr f)))
 195      (go loop)))
 196
 197 (defun cpbg0 (tr f)
 198   (prog (m f1 f2 g alc trm)
 199      (setq m 0)
 200      (cond ((and (not (numberp (caddr tr))) (alg (caddr tr)))
 201             (setq alc (painvmod (caddr tr)) tr (ptimes alc tr)))
 202            (t (setq alc 1.)))
 203      bk   (cond ((pcoefp f)
 204                  (return (cond ((and (null f1) (null f2))
 205                 ;; NOTE TO TRANSLATORS: MEANING OF NEXT MESSAGE IS OBSCURE
 206                                 (merror (intl:gettext "CPBG0: wrong trace.")))
 207                                (t
 208                                 (cons f1 f2)))))
 209                 ((equal (cadr f) 1)
 210                  (return (cons (cons f f1) f2)))
 211                 ((equal m modulus)
 212                  (return (cons f1 (cons f f2)))))
 213      (setq trm (pdifference tr (ptimes m alc)))
 214      (setq g (pgcdu trm f))
 215      (cond ((or (numberp g) (and *alpha* (alg g)))
 216             (incf m)
 217             (go bk)))
 218      (setq f (car (pmodquo f g)))
 219      (cond ((equal (cadr g) 1) (setq f1 (cons g f1)))
 220            (t (setq f2 (cons g f2))))
 221      (go bk)))
 222
 223 (defun cpol2p (p var)
 224   (prog((i 0) ans)
 225      (setq  p (nreverse p))
 226      loop (cond ((null p) (return (cons var ans)))
 227                 ((equal 0 (car p)) nil)
 228                 (t (setq ans (cons i (cons (car p) ans)))))
 229      (setq p (cdr p)
 230            i (1+ i))
 231      (go loop)))
 232
 233 (defun tracemod (v)
 234   (prog (ans tr qlarge term)
 235      (setq ans 0
 236            tr (nreverse trl*)
 237            trl* nil)
 238      (cond ((and (atom (caar tr)) (not (numberp (caar tr))))
 239             (setq qlarge t)))
 240      loop (when (null tr) (return ans))
 241      (setq term (if qlarge
 242                     (car tr)
 243                     (cpol2p (car tr) v))
 244            tr (cdr tr))
 245      (setq ans (pplus ans term))
 246      (setq trl* (cons term trl*))
 247      (go loop)))
 248
 249 (defun otracemod (term q m prime)
 250   (prog (ans i)
 251      (setq ans term
 252            i 1
 253            trl* (list term))
 254      loop (when (equal i m) (return ans))
 255      (setq ans (pplus ans (setq term (pexptmod term prime q))))
 256      (setq trl* (cons term trl*))
 257      (incf i)
 258      (go loop)))
 259
 260 (defun tracemod0 (q i)
 261   (prog (l ans a dl)
 262      (cond ((= i 0) (return (if trl*
 263                                 (tracemod (car q))
 264                                 (otracemod var q mm* modulus))))
 265            (trl* (setq dl trl*
 266                        trl* (mapcar #'(lambda(x)
 267                                         (cons (car x) (pgcd1 (cdr x) (cdr q)))) trl*))))
 268      (cond (tra* (go tag))
 269            (t (setq l (cdr trl*)
 270                     tra* (list *alpha*)
 271                     a *alpha*)))
 272      loop (when (null l) (go tag))
 273      (setq l (cdr l)
 274            a (pexpt a modulus)
 275            tra* (cons a tra*))
 276      (go loop)
 277      tag
 278      (setq ans (tracemod1 i tra* trl*))
 279      (when dl (setq trl* dl))
 280      (return ans)))
 281
 282 (defun tracemod1 (n a l)
 283   (prog (ans)
 284      (setq ans 0)
 285      loop (when (null l) (return ans))
 286      (setq ans (pplus ans (ptimes (pexpt (car a) n) (car l))))
 287      (setq l (cdr l)
 288            a (cdr a))
 289      (go loop)))
 290
 291 ;; The way arrays are manipulated has been changed to make this code reentrant.
 292 ;; Previously, arrays were kept on the array properties of symbols.  Now, the
 293 ;; symbols are bound to the arrays, so they can be rebound on re-entry.
 294 ;; The ANOTYPE, INVC, and FCTC arrays are set up in RAT;FACTOR.
 295
 296 (declare-top (special anotype invc fctc))
 297
 298 (defmacro a (row col)
 299   `(aref anotype ,row ,col))
 300
 301 (defmacro invc (ind)
 302   `(aref invc ,ind))
 303
 304 (defmacro fctc (ind)
 305   `(aref fctc ,ind))
 306
 307 (defun cptomexp (p m u n)
 308   (prog (b1 b2 j n-1 i l)
 309      (setq b1 (x**q1 (list (car u) 1 1) u m p))
 310      (cond ((equal (cdr b1) '(1 1))
 311             (setq split* t)
 312             (return nil)))
 313      (setq b2 b1 j 1 n-1 (1- n))
 314      (go tag1)
 315      tag (incf j)
 316      (when (= j n) (return nil))
 317      (setq b1 (pmodrem(ptimes b1 b2) u))
 318      tag1 (setq l (p2cpol b1 n-1) i n-1)
 319      sharp2   (when (null l) (go on))
 320      (setf (a j i) (car l))
 321      (setq l (cdr l))
 322      (setq i (1- i))
 323      (go sharp2)
 324      on   (setf (a j j) (pdifference (a j j) 1))
 325      (go tag)))
 326
 327 (defvar thr* 100)
 328
 329 (defun cptom (p m u n)
 330   (prog (( q (expt p m)) l s *xn (j 0) (i 0) ind n-1)
 331      (declare (special q i j))
 332      (setq  n-1 (1- n))
 333      (when (> q thr*) (return (cptomexp p m u n)))
 334      loop (incf j)
 335      (cond ((= j n) (return nil))
 336            (ind (go sa))
 337            (t
 338             (setq *xn (mapcar #'pminus (p2cpol (cddr u) n-1))
 339                   s (x**q (p2cpol(list var 1 1) n-1) p m)
 340                   ind t)))
 341      (go st)
 342      sa (cptimesxa s q)
 343      st (cond ((and (= j 1)
 344                     (equal '(1 0) (last s 2))
 345                     (= 1 (length (delete 0 (copy-tree s) :test #'equal))))
 346                (return (setq split* t))))
 347      (setq l s)
 348      (setq i n-1)
 349      sharp2   (when (null l) (go on))
 350      (setf (a j i) (car l))
 351      (setq l (cdr l))
 352      (decf i)
 353      (go sharp2)
 354      on   (setf (a j j) (pdifference (a j j) 1))
 355      (go loop)))
 356
 357 (defun cptimesxa (p i)
 358   (prog (xn q lc)
 359    ag    (when (= i 0) (return p))
 360    (setq xn *xn
 361          q p
 362          lc (car p))
 363    loop (cond ((cdr q)
 364                (rplaca q (pplus (cadr q) (ptimes lc (car xn))))
 365                (setq q (cdr q) xn (cdr xn)))
 366               (t (rplaca q (ptimes lc (car xn)))
 367                  (decf i)
 368                  (go ag)))
 369    (go loop)))
 370
 371 (defun x**q (x p m)
 372   (prog ((i 1)  (pp 1) (d 0))
 373      (setq i 1 trl* (list x) pp 1)
 374      loop (when (= i m) (return (cptimesxa x (- (* pp p) pp))))
 375      (setq d pp)
 376      (cptimesxa x (- (setq pp (* pp p)) d))
 377      (setq trl* (cons (copy-tree x) trl*))
 378      (incf i)
 379      (go loop)))
 380
 381 (defun cmnull (n)
 382   (prog (nullsp (sub1n (1- n)) mone (k 1) (j 0) (s 0) nullv (i 0) vj m aks)
 383      (setq mone (cmod -1))
 384      sharp    (when (> i sub1n) (go on))
 385      (setf (fctc i) -1)
 386      (setf (invc i) -1)
 387      (incf i)
 388      (go sharp)
 389      on   (setq k 1 nullsp (list 1))
 390      n2   (when (> k sub1n) (return nullsp))
 391      (setq j 0)
 392      n3a  (cond ((> j sub1n) (go null))
 393                 ((or (equal (a k j) 0) (> (fctc j) -1))
 394                  (incf j)
 395                  (go n3a)))
 396      (setf (invc k) j)
 397      (setf (fctc j) k)
 398      (setq m (a k j))
 399      (setq m (ptimes mone (painvmod m)))
 400      (setq s k)
 401      sharp1   (when (> s sub1n) (go on1))
 402      (setf (a s j) (ptimes m (a s j)))
 403      (incf s)
 404      (go sharp1)
 405      on1
 406      (setq s 0)
 407      sharp2   (when (> s sub1n) (go nextk))
 408      (cond ((= s j) nil)
 409            (t (prog (i)
 410                  (setq aks (a k s))
 411                  (setq i k)
 412                  sharp3   (when (> i sub1n) (return nil))
 413                  (setf (a i s) (pplus (a i s) (ptimes (a i j) aks)))
 414                  (incf i)
 415                  (go sharp3))))
 416      (incf s)
 417      (go sharp2)
 418      null (setq nullv nil)
 419      (setq s 0)
 420      sharp4   (cond ((> s sub1n) (go on4))
 421                     ((= s k) (setq nullv (cons s (cons 1 nullv))))
 422                     ((> (invc s) -1)
 423                      (setq vj (a k (invc s)))
 424                      (cond ((equal vj 0) nil)
 425                            (t (setq nullv (cons s (cons vj nullv)))))))
 426      (incf s)
 427      (go sharp4)
 428      on4  (cond ((equal (car nullv) 0) (setq nullv (cadr nullv)))
 429                 ((setq nullv (cons var nullv))))
 430      (setq nullsp (cons nullv nullsp))
 431      nextk (incf k)
 432      (go n2)))