src/nalgfa.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 1980 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module nalgfa)
  14
  15 (declare-top (special vlist *nosplitf *algvar *denom *num *ans
  16                       algfac* $nalgfac alpha))
  17
  18 (load-macsyma-macros rzmac ratmac)
  19
  20 (defun new-alg ()
  21   (newsym ($gensym "alg")))
  22
  23
  24 (defun psqfrp (p var)
  25   (zerop (pdegree (pgcd p (pderivative p var)) var)))
  26
  27 (defun pgsubst (val var p) ;;generalized psubst substitutes any
  28   (cond ((pcoefp p) p) ;;expression for any var in p
  29         ((eq var (car p))
  30          (cond ((pzerop val)
  31                 (ptterm (cdr p) 0))
  32                ((do ((ld (cadr p) (car a))
  33                      (a (cdddr p) (cddr a))
  34                      (ans (caddr p)
  35                           (pplus
  36                            (ptimes ans (pexpt val (- ld (car a))))
  37                            (cadr a))))
  38                     ((null a) (ptimes ans (pexpt val ld)))))))
  39         ((pointergp var (car p)) p)
  40         ((do ((a (cdddr p) (cddr a))
  41               (ans (ptimes (list (car p) (cadr p) 1)
  42                            (pgsubst val var (caddr p)))
  43                    (pplus ans
  44                           (ptimes (list (car p) (car a) 1)
  45                                   (pgsubst val var (cadr a))))))
  46              ((null a) ans)))))
  47
  48 (defun pvsubst (nvar ovar p)
  49   (cond ((or (pcoefp p) (pointergp ovar (car p))) p)
  50         ((eq ovar (car p))
  51          (cons nvar (cdr p)))
  52         (t (pgsubst (make-poly nvar) ovar p))))
  53
  54 (defun ordervar (var l)
  55   (let ((mvar (lmainvar l)))
  56     (cond ((null mvar) l)
  57           ((null (pointergp mvar var)) (cons var l))
  58           ((let ((newvar (gensym)))
  59              (setq genvar (append genvar (list newvar)))
  60              (setf (symbol-plist newvar) (symbol-plist var))
  61              (setf (symbol-value newvar) (1+ (symbol-value mvar)))
  62              (cons newvar (mapcar #'(lambda (p) (pvsubst newvar var p)) l)))))))
  63
  64 (defun lmainvar (l) ;;main var of list of poly's
  65   (do ((l l (cdr l))
  66        (v))
  67       ((null l) v)
  68     (cond ((pcoefp (car l)))
  69           ((null v) (setq v (caar l)))
  70           ((pointergp (caar l) v)
  71            (setq v (caar l))))))
  72
  73 (defun presult (p1 p2 var) ;;change call in algsys?
  74   (let ((genvar genvar))
  75     (setq var (ordervar var (list p1 p2))
  76           p1 (cadr var)
  77           p2 (caddr var)
  78           var (car var))
  79     (cond ((zerop (pdegree p1 var))
  80            (cond ((zerop (pdegree p2 var)) 1)
  81                  ((pexpt p1 (cadr p2)))))
  82           ((zerop (pdegree p2 var))
  83            (pexpt p2 (cadr p1)))
  84           ((resultant p1 p2)))))
  85
  86 (defun pcoefvec (p)
  87   (cond ((pcoefp p) (list p))
  88         ((do ((l)
  89               (i (cadr p) (1- i))
  90               (p (cdr p)))
  91              ((signp l i) (nreverse l))
  92            (push (cond ((and p (= (car p) i))
  93                         (prog1 (cadr p) (setq p (cddr p))))
  94                        ( 0 ))
  95                  l)))))
  96
  97 (defun algtrace1 (bvec tvec)
  98   (do ((i (- (length tvec) (length bvec)) (1- i)))
  99       ((zerop i) (algtrace* bvec tvec))
 100     (setq bvec (cons 0 bvec))))
 101
 102 (defun algtrace* (bvec tvec)
 103   (do ((b bvec (cdr b))
 104        (tr (car (last bvec))
 105            (pplus tr (car (last b)))))
 106       ((null (cdr b)) tr)
 107     (or (pzerop (car b))
 108         (do ((l (cdr b) (cdr l))
 109              (tv tvec (cdr tv)))
 110             ((null l))
 111           (rplaca l (pdifference (car l) (ptimes (car b) (car tv))))))))
 112
 113 (defun algtrace (r p)
 114   (cond ((eq (caar r) (car p))
 115          (ratreduce (algtrace1 (pcoefvec (car r))
 116                                (cdr (pcoefvec p)))
 117                     (cdr r)))
 118         ((ratreduce (pctimes (cadr p) (car r))
 119                     (cdr r)))))
 120
 121 (defmfun $algtrace (r p var)
 122   (let ((varlist (list var))
 123         (genvar nil))
 124     (rdis* (algtrace (rform r) (car (rform p))))))
 125
 126 (defun good-form (l) ;;bad -> good
 127   (do ((l l (cddr l))
 128        (ans))
 129       ((null l) (nreverse ans))
 130     (push (cons (cadr l) (car l)) ans)))
 131
 132 (defun bad-form (l) ;;good -> bad
 133   (mapcar #'(lambda (q) (list (cdr q) (car q))) l))
 134
 135 (defmfun $algfac (a1 &optional (a2 nil a2?) (a3 nil a3?))
 136   (if a3?
 137       ($pfactoralg a1 a2 a3)
 138       (let ((varlist))
 139         (cond (a2?
 140                (newvar a2)
 141                (if (alike1 a2 (car varlist))
 142                    ($pfactoralg a1 nil a2)
 143                    ($pfactoralg a1 a2 (car (last varlist)))))
 144               (t
 145                (newvar a1)
 146                  (setq varlist (mapcan #'(lambda (q) (if (algpget q) (list q) nil)) varlist))
 147                  (cond ((= (length varlist) 1)
 148                         ($pfactoralg a1 nil (car varlist)))
 149                        ((> (length varlist) 1)
 150                         ;; MEANING OF NEXT MESSAGE IS UNCLEAR
 151                         (merror (intl:gettext "algfac: too many algebraics.")))
 152                        (t
 153                         ;; MEANING OF NEXT MESSAGE IS UNCLEAR
 154                         (merror (intl:gettext "algfac: no algebraics.")))))))))
 155
 156 (defmfun $pfactoralg (f p alg)
 157   (let ((varlist (list alg))
 158         (genvar) (vlist) (tellratlist) ($ratfac)
 159         ($gcd '$algebraic)
 160         ($algebraic) ($ratalgdenom t)
 161         (*denom 1) (*num 1) (*ans))
 162     (cond ((and (null p) (radfunp alg t)) (newvar (cadr alg)))
 163           (t (newvar p)))
 164     (newvar1 f)
 165     (cond ((null vlist) (merror (intl:gettext "pfactoralg: attempt to factor a constant."))))
 166     (setq varlist (nconc varlist (sortgreat vlist)))
 167     (cond (p (setq p (cadr (ratrep* p)))
 168              (push (cons alg (mapcar #'pdis (cdr p)))
 169                    tellratlist))
 170           (t (setq p (algpget alg))
 171              (setq p (pdifference
 172                       (pexpt (cadr (ratrep* alg)) (car p))
 173                       (cadr p)))))
 174     (setq $algebraic t)
 175     (setq f (cadr (ratrep* f)))
 176     (setq f (pfactoralg1 f p 0))
 177     (cons '(mtimes)
 178           (cons (rdis (ratreduce *num *denom))
 179                 (mapcar 'pdis f)))))
 180
 181 (declare-top (special adn*)) ;also redefine fact5 to call nalgfac correctly
 182
 183 (defun nalgfac (p mp)
 184   (let ((*num 1) (*denom 1) (*ans) (algfac*) ($nalgfac)
 185         ($algebraic t) ($gcd '$algebraic))
 186     (setq p (pfactoralg1 p mp 0))
 187     (setq adn* (* adn* *denom))
 188     (cond ((equal *num 1) p)
 189           (t (cons *num p)))))
 190
 191 (setq *nosplitf t)
 192
 193 (defun pfactoralg1 (f p ind)
 194   (cond ((pcoefp f) (setq *num (ptimeschk f *num)) *ans)
 195         ((= (cadr f) 1) (setq f (pshift f (car p) ind))
 196          (push (algnormal f) *ans)
 197          (setq f (rquotient f (car *ans))
 198                *denom (ptimeschk (cdr f) *denom)
 199                *num (ptimeschk (car f) *num))
 200          *ans)
 201         ((equal (cdr f) (cdr p))
 202          (push (pdifference (make-poly (car f)) (make-poly (car p))) *ans)
 203          (setq f (rquotient f (car *ans))
 204                *denom (ptimeschk (cdr f) *denom))
 205          (pfactoralg1 (car f) p ind))
 206         ((zerop (pdegree f (car p)))
 207          (mapc #'(lambda (q)
 208                    (if (pcoefp q) nil
 209                      (pfactoralg1 (pshift q (car p) -1) p (1+ ind))))
 210                (let (($algebraic nil)
 211                      ($gcd (car *gcdl*)))
 212                  (pfactor1 f)))
 213          *ans)
 214         (t (do ((l (let (($algebraic nil)
 215                          ($gcd (car *gcdl*)))
 216                      (pfactor (algnorm f p)))
 217                    (cddr l))
 218                 (polys)
 219                 (temp)
 220                 (alg (car p)))
 221                ((null l)
 222                 (setq *num (ptimeschk f *num))
 223                 (mapc #'(lambda (q) (pfactoralg1
 224                                      (pshift q alg -1) p (1+ ind)))
 225                       polys)
 226                 *ans)
 227                (cond ((pcoefp (car l)) nil)
 228                      (t (setq temp (cond ((null (cddr l)) f)
 229                                          (t (pgcd f (car l)))))
 230                         (cond ((pcoefp temp) nil)
 231                               ((= (cadr temp) 1)
 232                                (setq temp (algnormal temp))
 233                                (push (pshift temp alg ind) *ans))
 234                               ((= (cadr l) 1)
 235                                (setq temp (algnormal temp))
 236                                (push (pshift temp alg ind) *ans)
 237                                (or *nosplitf
 238                                    (setq *nosplitf
 239                                          (list (car l) temp ind))))
 240                               (t (push temp polys)))
 241                         (setq f (rquotient f temp)
 242                               *denom (ptimeschk (cdr f) *denom)
 243                               f (car f)))) ))))
 244
 245 (defun pshift (f alg c)
 246   (if (= c 0) f
 247       (pgsubst (pplus (make-poly (car f)) (pctimes c (make-poly alg)))
 248                (car f) f)))
 249
 250
 251
 252 (defmfun $splitfield (p var)
 253   (let ((varlist)
 254         (genvar)
 255         (genpairs)
 256         (*algvar)
 257         ($gcd '$algebraic))
 258     (newsym var)
 259     (setq *algvar (caar (new-alg)))
 260     (setq p (psplit-field (cadr (ratf p))))
 261     (cons
 262      '(mlist)
 263      (cons (pdis* (car p))
 264            (mapcar 'rdis* (cdr p))))))
 265
 266 (defun psplit-field (p)                 ;modresult?
 267   (let ((l (mapcar #'(lambda (q) (psplit-field1 (cdr q)))
 268                    (good-form (pfactor p)))) ;don't normalize lcfs?
 269         ($algebraic t))
 270     (if (null (cdr l)) (car l)
 271         (do ((l l (cdr l))
 272              (prim) (zeroes) (temp))
 273             ((null l) (cons (or prim '|$splits in q|) zeroes))
 274           (cond ((eq (caar l) 'linear)
 275                  (setq zeroes (cons (cdar l) zeroes)))
 276                 ((null prim)
 277                  (setq prim (caar l)
 278                        zeroes (nconc (cdar l) zeroes)))
 279                 ((setq temp
 280                        (primelmt (cons (car p) (cdr prim))
 281                                  (cons (car p) (cdaar l))
 282                                  *algvar)
 283                        zeroes
 284                        (nconc
 285                         (mapcar
 286                          #'(lambda (q)
 287                              (ratgsubst (cadddr temp) (caaar l) q))
 288                          (cdar l))
 289                         (mapcar
 290                          #'(lambda (q)
 291                              (ratgsubst (caddr temp) (car prim) q))
 292                          zeroes))
 293                        prim (car temp))))))))
 294
 295 (defun plsolve (p)
 296   (ratreduce (ptterm (cdr p) 0)
 297              (pminus (ptterm (cdr p) 1))))
 298
 299
 300 (defun psplit-field1 (p)
 301   ;;returns minimal poly and list of roots
 302   ;;p must be square free
 303   (*bind* ((minpoly (cons *algvar (cdr p)))
 304            (zeroes) ($algebraic t)
 305            ($ratalgdenom t))
 306           (if (equal (cadr p) 1) (return (cons 'linear (plsolve p))))
 307           (do ((polys (list p) )
 308                (nminpoly)
 309                (*nosplitf nil nil)
 310                (alpha (cons (make-poly (car minpoly)) 1)))
 311               ((null polys)
 312                (cons minpoly zeroes))
 313             (push alpha zeroes)
 314             (putprop (car minpoly) (cdr minpoly) 'tellrat)
 315             (rplaca polys
 316                     (car
 317                      (rquotient (pctimes (cdr alpha) (car polys))
 318                                 (pdifference
 319                                  (pctimes (cdr alpha) (pget (caar polys)))
 320                                  (car alpha)))))
 321             (setq polys
 322                   (mapcan
 323                    #'(lambda (q)
 324                        (cond ((equal (cadr q) 1) ;;linear factor
 325                               (push (plsolve q) zeroes)
 326                               nil) ;;flush linear factor
 327                              ((list q))))
 328                    (mapcan #'(lambda (q)
 329                                (let ((*ans) (*num 1) (*denom 1))
 330                                  (nreverse (pfactoralg1 q minpoly 0))))
 331                            polys)))
 332             (when *nosplitf
 333               (setq nminpoly (car *nosplitf)
 334                     *nosplitf (cdr *nosplitf))
 335               (putprop *algvar (cdr nminpoly) 'tellrat)
 336               (let ((beta
 337                      (plsolve (pgcd (cons (caar *nosplitf) (cdr minpoly))
 338                                     (exchangevar (car *nosplitf) *algvar)))))
 339                 (setq alpha (ratplus (cons (make-poly *algvar) 1)
 340                                      (rattimes (cons (- (cadr *nosplitf)) 1)
 341                                                beta t)))
 342                 (setq zeroes
 343                       (mapcar
 344                        #'(lambda (q) (ratgsubst beta (car minpoly) q))
 345                        zeroes))
 346                 (setq polys
 347                       (mapcar
 348                        #'(lambda (q) (car (rgsubst beta (car minpoly) q)))
 349                        polys))
 350                 (setq minpoly
 351                       (cons *algvar (cdr nminpoly))))))))
 352
 353 (defun exchangevar (poly var)
 354   (let ((newvar (gensym))
 355         (ovar (car poly)))
 356     (setf (symbol-value newvar) (1+ (eval ovar)))
 357     (pvsubst ovar newvar
 358              (pvsubst var ovar
 359                       (pvsubst newvar var poly)))))
 360
 361 (defun rgsubst (val var p) ;;generalized psubst substitutes any
 362   (cond ((pcoefp p)
 363          (cons p 1)) ;;expression for any var in p
 364         ((eq var (car p))
 365          (cond ((pzerop val)
 366                 (cons (ptterm (cdr p) 0) 1))
 367                ((do ((ld (cadr p) (car a))
 368                      (a (cdddr p) (cddr a))
 369                      (ans (cons (caddr p) 1)
 370                           (ratplus
 371                            (rattimes ans
 372                                      (ratexpt val
 373                                               (- ld (car a)))
 374                                      t)
 375                            (cons (cadr a) 1))))
 376                     ((null a) (rattimes ans (ratexpt val ld) t))))))
 377         ((pointergp var (car p)) (cons p 1))
 378         (t (let ((newsym (gensym)))
 379              (setf (symbol-value newsym) (1+ (symbol-value (car p))))
 380              (rgsubst val newsym (pvsubst newsym var p))))))
 381
 382 (defun ratgsubst (val var rat)
 383   (ratquotient (rgsubst val var (car rat))
 384                (rgsubst val var (cdr rat))))
 385
 386 (defun algnorm (f p)
 387   (presult f p (car p)))
 388
 389 (defmfun $algnorm (r p var)
 390   (let ((varlist (list var))
 391         (genvar))
 392     (setq r (ratf r)
 393           p (cadr (ratf p)))
 394     (rdis* (cons (algnorm (cadr r) p)
 395                  (algnorm (cddr r) p)))))
 396
 397 (defun sqfrnorm (f p fvar) ;;f must be sqfr, p is minpoly, fvar # pvar
 398   (*bind* ((pvar (car p)))
 399           (setq f (cdr (ordervar pvar (list f p))) ;;new main var will be car of p
 400                 p (cadr f) f (car f))  ;make mainvar of f = mainvar(p)
 401           (do ((i 0 (1+ i))
 402                (dif (pdifference (make-poly fvar) (make-poly (car p))))
 403                (f f (pgsubst dif fvar f))
 404                (res))
 405               ((and (eq (car f) (car p))
 406                     (setq res (primpart (algnorm f p)))
 407                     (psqfrp res fvar))
 408                (list res
 409                      (*bind* (($algebraic t))           ;;;modified f
 410                              (putprop pvar (cdr p) 'tellrat)
 411                              (pvsubst pvar (car p) f))
 412                      (car p)
 413                      p
 414                      i)))))
 415
 416 (defun primelmt (a b gvar &aux ($algebraic nil))
 417   ;;a is a poly with coeff's in k(b)
 418   ;;gvar is new variable
 419   (let ((norm (sqfrnorm (cons gvar (cdr a)) b gvar))
 420         (alpha) (beta) ($ratalgdenom t))
 421     (rplaca norm (primpart (car norm)))
 422     (putprop gvar (cdar norm) 'tellrat)
 423     (setq $algebraic t
 424           beta (pgcd (cadddr norm)
 425                      (pvsubst (caddr norm)
 426                               (car b)
 427                               (cadr norm))))
 428     (setq beta (plsolve beta)
 429           alpha (ratplus (cons (make-poly gvar) 1)
 430                          (rattimes (cons (- (cadddr (cdr norm))) 1)
 431                                    beta t)))
 432     (list (car norm) ;;minimal poly
 433           (pplus (make-poly a) ;;new prim elm in old guys
 434                  (list (car b) 1 (- (cadddr (cdr norm)))))
 435           alpha beta))) ;;in terms of gamma
 436
 437 ;; discriminant of a basis
 438
 439 (defmfun $bdiscr (&rest args)
 440   (let ((varlist) (genvar))
 441     (xcons (bdiscr (mapcar #'rform (butlast args))
 442                    (car (rform (car (last args)))))
 443            (list 'mrat 'simp varlist genvar))))
 444
 445 (defun bdiscr (l minp)
 446   (det (mapcar #'(lambda (q1)
 447                    (mapcar #'(lambda (q2) (algtrace (ptimes q1 q2) minp)) l))
 448                l)))